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Abstracts 



English abstract 

In interacting galaxies, strong tidal forces disturb the global morphology of the pro- 
genitors and give birth to the long stellar, gaseous and dusty tails often observed. In 
addition to this destructive effect, tidal forces can morph into a transient, protective 
setting called compressive mode. Such modes then shelter the matter in their midst by 
increasing its gravitational binding energy. 

This thesis focuses on the study of this poorly known regime by quantifying its prop- 
erties thanks to numerical and analytical tools applied to a spectacular merging system 
of two galaxies, commonly known as the Antennae galaxies. iV-body simulations of this 
pair yield compressive modes in the regions where observations reveal a burst of star 
formation. Furthermore, characteristic time- and energy scales of these modes match 
well those of self-gravitating substructures such as star clusters and tidal dwarf galaxies. 
Comparisons with star formation rates derived from hydrodynamical runs confirm the 
correlation between the location of compressive modes and sites where star formation is 
likely to show enhanced activity. Altogether, these results suggest that the compressive 
modes of tidal fields plays an important role in the formation and evolution of young 
clusters, at least in a statistical sense, over a lapse of ~ 10 million years. Preliminary 
results from simulations of stellar associations highlight the importance of embedding the 
clusters in the evolving background galaxies to account precisely for their morphology 
and internal evolution. 

These conclusions have been extended to numerous configurations of interacting galax- 
ies and remain robust to a variation of the main parameters that characterize a merger. 
We report however a clear anti-correlation between the importance of the compressive 
mode and the distance between the galaxies. Further studies including hydrodynamics 
are now underway and will help pin down the exact role of the compressive mode on the 
formation and later survival of star clusters. Early comparisons with such computations 
suggest that compressive modes act as catalysts or triggers of star formation. 
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Resume en francais 

Dans les galaxies en interaction, de colossales forces de maree perturbent la morpholo- 
gie des progeniteurs pour engendrer les longs bras d'etoiles, gaz et poussieres que Ton 
observe parfois. En plus de leur effet destructeur, les forces de maree peuvent, dans 
certain cas, se placer dans une configuration protectrice appelee mode compressif. De 
tels modes protegent alors la matiere en leur sein, en augmentant son energie de liaison. 

Cette these se concentre sur l'etude de ce regime peu connu en quantifiant ses pro- 
prietes grace a des outils numeriques et analytiques appliques a un spectaculaire systeme 
de galaxies en fusion, communement appele les Antennes. Des simulations iV-corps de 
cette paire de galaxies montrent la presence de modes compressifs dans les regions ou les 
observations revelent un sursaut de formation stellaire. De plus, les temps et energies 
caracteristiques de ces modes correspondent a ceux de la formation de sous-structures 
autogravitantes telles que des amas stellaires et des naines de maree. Des comparaisons 
avec les taux de formation stellaire derives de simulations hydro dynamiques confirment la 
correlation entre les positions des modes compressifs et les sites ou la formation des etoiles 
est certainement amplifiee. Mis bout-a-bout, ces resultats suggerent que les modes com- 
pressifs des champs de maree jouent un role important dans la formation et revolution 
des jeunes amas, au moins d'un point de vue statistique, sur une echelle de temps de 
l'ordre de dix millions d'annees. Des resultats preliminaries de simulations d'associations 
stellaires soulignent l'importance de plonger les amas dans leur environnement galactique 
en evolution, pour tenir compte precisement de leur morphologie et evolution interne. 

Ces conclusions ont ete etendues a de nombreuses configurations d'interaction et 
restent robustes aux variations des principaux parametres caracterisant les paires de 
galaxies. Nous notons cependant une nette anti-correlation entre l'importance du mode 
compressif et la distance entre ces galaxies. De nouvelles etudes incluant les aspects 
hydrodynamiques sont maintenant en cours et aideront a preciser le role exact du mode 
compressif dans la formation et la survie des amas d'etoiles. Les premieres compara- 
isons avec de telles simulations suggerent que les modes compressifs agissent en tant que 
catalyseurs ou amorces de la formation stellaire. 
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Deutsch Zusammenfassung 

Starke Gezeitenkrafte in wechselwirkenden Galaxien storen die Morphologie dieser 
Systeme und fordern die Entstehung ausgedehnter und haufig beobachteter Filamente 
aus Sternen, Gas und Staub. Neben diesem zerstorerischen Effekt konnen Gezeitenkrafte 
auch zu einem voriibergehenden stabilisierenden Zustand, den sogenannten kompressiven 
Moden fiihren. Durch die Erhohung der gravitativen Bindungsenergie der vorhandenen 
Materie wird diese dann von aueren gravitativen Einflassen abgeschirmt. 

Die vorliegende Arbeit beschaftigt sich mit diesen wenig untersuchten Zustanden 
durch Quantifizierung ihrer Eigenschaften mittels numerischer und analytischer Metho- 
den, angewandt auf ein spektakulares System verschmelzender Galaxien, bekannt als die 
Antennengalaxien. N-K6rper Simulationen dieses Galaxienpaares ergeben kompressive 
Moden in denselben Regionen wo Beobachtungen eine erhohte Sternentstehung zeigen. 
Die charakteristischen Zeit- und Energieskalen dieser Moden ahneln auBerdem stark de- 
nen selbstgravitierender Substrukturen wie Sternhaufen oder sogenannten tidal dwarfs. 
Vergleiche mit Sternentstehungsraten aus hydrodynamischen Simulationen bestatigen die 
Korrelation zwischen der Lage der kompressiven Moden und Bereichen mit erhohter Ster- 
nentstehung. Zusammenfassend ist zu sagen, dass diese Resultate darauf hinweisen, dass 
kompressive Moden von Gezeitenfeldern statistisch betrachtet eine wichtige Rolle bei der 
Entstehung und Entwicklung junger Sternhaufen in einem Zeitraum von 10 Millionen 
Jahren spielen. Vorlaufige Resultate von Simulationen stellarer Assoziationen zeigen die 
Bedeutung der Einbettung dieser Haufen in die sich entwickelnden Muttergalaxien, um 
deren Morphologie und interne Entwicklung zu begrunden. 

Diese Schlussfolgerungen wurden auf zahlreiche Konfigurationen wechselwirkender 
Galaxien erweitert und liefern bei Variation der charakteristischen Parameter ver- 
schmelzender Galaxien das gleiche Ergebnis. Es ist jedoch eine klare Anti-Korrelation 
zwischen der Wichtigkeit der kompressiven Moden und dem Abstand der Galaxien 
zueinander zu erkennen. Weitere hydro dynamische Studien sind nun in vollem Gange und 
werden dazu beitragen, den genauen Einfluss der kompressiven Moden auf die Entste- 
hung und dem spaterem Uberleben von Sternhaufen festzulegen. Aktuelle Vergleiche 
mit solchen Berechnungen weisen darauf hin, dass kompressive Moden als Katalysatoren 
oder Ausloser von Sternentstehung angesehen werden konnen. 
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Preamble 



Prologue 

One of the particularities of this PhD is that it has been co-supervised. Indeed, I 
worked with two different advisors in two different institutes in two different countries. 
Advantages of such an organization are indubitably the richness of the ideas and advices 
I collected from both sides. Of course, sometimes I had to deal with points of view that 
may be incompatible, but generally, this helped me to learn how to make decisions and 
to stand up for them. 




Left: Das Institut fur Astronomie der Universitat Wien (Osterreich). 
Right: L'Observatoire Astronomique de Strasbourg (France). 



This also requires much more financial means than a classical PhD and a good flex- 
ibility in the schedule. Fortunately for me, both Institutes and Universities supported 
this collaboration and provided the resources to make it a wonderful experience. 

The work presented here begun within the framework of a previous PhD thesis "Ap- 
proche numerique de la dynamique et de revolution stellaires appliquees a la fusion 
galactique" 3 , led by Jean-Julien Fleck under the supervision of Christian Boily at the 
Observatoire de Strasbourg (2004-2007). Their pioneer methods have been optimized to 
extend their conclusions about the tidal field and the star clusters. 



a "Numerical approach of stellar dynamic and evolution applied to galactic mergers" . The manuscript 
(in french) can be downloaded: http://tel.archives-ouvertes.fr/tel-00184822 [S 1 
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How to read this document 



This manuscript aims to present the scientific motivations, methods, results, questions 
answered and questions raised during my PhD. It has been conceived to introduce this 
work but also to archive the results and the mathematical derivations. That is, the 
General Introduction gives a broad review about the major topics developed in the next 
Chapters, while the Appendices gather technical and analytical work, and give some 
hints for the derivation of several equations. 



Box 1 

Supplementary "boxes" like this one give more details on technical and/or 
mathematical points. They can be skipped during the first reading. 



In the electronic version of this manuscript, hyperlinks point to all referenced figures, 
tables, equations and sections. Within a citation, clicking on the year will direct you to 
the General Bibliography where you will find a link to the ADS page (which will open 
in you browser) of the publication. In addition, it is possible to click on the acronyms, 
and go where they have been defined for the first time. In this (illustrative) example: 

Once you have seen Figure 11.12 (page 73), you may be 
interested in looking at Table II. 1 and find out more 
about the SFE in Section III. 2.2 or in Renaud et al. 
(2009). 

you can click on "Figure 11.12", "Table II.l", "SFE", "Section III.2.2" and "2009". Don't 
forget to use the Back function of your favorite PDF reader to go back to the link! 



The author 

somewhere in the Alps, between Vienna and Strasbourg 
Wednesday, 17 February 2010 (1:15 am) 
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General Introduction 



Overview 

This Chapter introduces briefly the historical ideas of 
galaxy mergers and star clusters, and the physical con- 
cepts we will refer to in the next parts. 

J 



1 Galaxy mergers 

1.1 How to classify the unclassified? 

For a long time after their discovery, galaxies were thought to be clouds or nebulae, 
named island universes, suggesting that they were isolated and without any interaction 
(see e.g. Hubble 1929). From Edwin Hubble (1936) to the recent Galaxy Zoo project 
(Raddick et al. 2007), an important body of work has aimed at classifying the galaxies 
according to their morphology, on the sky. Hubble's classification (or Hubble's pitch fork 
diagram, see Figure 1) distinguishes between two major families: the ellipticals and the 
spirals. 

This nomenclature has been very efficient for most of the galaxies and is still used 
nowadays. However, as the observational techniques and accuracy increased with time, 
more and more galaxies did not fit into this classification. Many peculiar morphologies 
showing extended arms lead to uncertainties, but it appears that a lot of them could be 
multiple galaxies, so closely packed that they would affect each other's shape. 

It rapidly became clear that the "island universes" were far from being isolated, as they 
are often associated in groups and clusters. Indeed, the environment plays a crucial role 
on the morphology of these galaxies. As two (or more) of them move close to each other, 
the gravitation brings them even closer and, if their relative velocity is low enough (i.e. 
less than the escape velocity), they finally merge into one single galaxy. Observations of 
various pairs of interacting galaxies suggested an evolutionary sequence (often called the 
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Figure 1: Hubble's classification of galaxies associates one or two letters for each family 
(E for the elliptical, S for the spirals, SB for the barred spirals) and, either a number to 
quantify either their axis ratio (ellipticals) or a letter for the opening angle of their spiral 
arms. (Picture adapted from Hubble Heritage Project website.) 



Toomre sequence; see Figure 2), in which each step represents a well defined dynamical 
stage (Toomre 1977; Joseph & Wright 1985; Hibbard & van Gorkom 1996; Laine et al. 
2003). 

Many phenomena are involved in such a process and thus mergers are wonderful 
laboratories for the astrophysicist. Among these phenomena, one counts tides, dynamical 
friction, shocks, star formation... An outstanding goal of this thesis is to weave them 
together to paint a more accurate picture of mergers. 



1.2 A telescope meets a computer 

Obviously, the dynamics of such systems is of prime importance to understand the 
history of these galaxies. Our knowledge of the motion of the progenitor galaxies and 
their initial morphology has to be as accurate as possible. On the one hand, the only 
information we dispose of is the integrated light of the galaxies, as viewed form the Earth. 
Therefore, we miss the third dimension. On the other hand, the radial velocity can be 
measured, and gives us an idea about the motion along the line of sight, but not along 
the other directions. Furthermore, even if astronomy has been one of the first sciences in 
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the history of mankind, accurate data has been collected over a very short period (few 
centuries) in regard to the timescales of galaxy mergers. And thus, we also miss the 
fourth dimension: time. 

Therefore the challenge of understanding interacting galaxies consists in understand- 
ing phenomena that have occurred far away, over a very long time, just by getting 
snapshots of the objects, from our perspective, at the present time. The role of observers 
is to provide the theoreticians with as rich data as possible, through photometry and 
spectroscopy. Then, we build models that (partially) reproduces the physics and try to 
interpret the observations. Whereas the models will never be perfect, they may allow to 
explore (up to a certain level) the unreachable quantities one needs to fully understand 
an object. 

The first step would be an analytical job. By making strong assumptions, one can 
mathematically derive some values of interest for simple models. However, the complexity 




Figure 2: The Tommre sequence (Toomre 1977) represents the supposed evolution of 
interacting galaxies. It starts with the early phases, when progenitors have just begun 
to interact, shows intermediate stages and finishes with the merger phase. From left to 
right, top: NGC 4038/39 (the Antennae), NGC 4676 (the Mice), NGC 3509, NGC 520; 
bottom: NGC 2623, NGC 3256, NGC 3921, NGC 7252. (Images from the Arp Atlas, 
montage by J. Hibbard.) 
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Figure 3: A light bulb and a photocell used by Holmberg as source and measurement 
tool of his proxy for gravitation. (Figure 1 of Holmberg 1941.) 



of the problem often forces the theoretician to make unrealistic hypotheses to simplify 
the mathematical developments. In the era of computers, we have a new way to solve the 
problems. The laws of physics that have been tested in laboratories are implemented as 
numerical treatments and applied to a set of initial conditions. The solution one obtains 
is no longer analytical, but numerical. At this point, the complexity becomes a second 
order problem while we have to face the precision versus speed issue. Indeed, once the 
physical laws have been implemented on the computer, the accuracy of the results is 
directly linked with the precision and the resolution of the computation. Usually, both 
can be increased at the cost of a longer calculation. Chapter I describes the numerical 
techniques used during this thesis and details the tests performed to validate the results. 

With telescopes providing data of higher and higher resolution, the theoreticians 
must adapt the techniques they implement and use more and more state-of-the-art (su- 
per) computers to recover the physics of smaller and smaller systems. 



1.3 When stars were light-bulbs 

In the case of galaxy mergers, simplifying assumptions are difficult to make: the 
distribution of matter offers no plane of symmetry, the masses are not negligible and 
the Keplerian laws of motion do not apply. All the usual simplifications are not possible 
because of the complexity level of the problem. Therefore, the use of simulations rapidly 
appeared vital to understand these systems. 

For his first attempt, Holmberg (1941) did not have computers. But he knew that both 
electromagnetic and gravitational interactions scale with the distance as r~ 2 . So he used 
37 light-bulbs to mimic gravitational sources and some galvanometers (see Figure 3) to 
measure the total amount of light at various points, i.e. the proxy for the gravitational 
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force. Therefore, he was able to determine the direction and intensity of his pseudo 
gravitational force and then move the sources according to the equations of motion. 

Two decades later, the first digital computer simulations used 100 particles to re- 
produce stellar systems in equilibrium to validate the integrators (von Hoerner 1960; 
Pfleiderer & Siedentopf 1961; Aarseth 1963; Eneev, Kozlov, & Sunyaev 1973) but the 
major breakthrough occurred in the seventies with the work of Toomre & Toomre (1972). 
They set disks of particles to reproduce two disk galaxies. Each galaxy consists in a point 
mass surrounded by test (i.e. mass-less) particles. Therefore, only the two-body prob- 
lem is to be solved, which can be done analytically through Kepler's laws of motion. 
This approach, called restricted simulation, allows to consider the effect of the potential 
(created by two bodies only) on a large number of test particles that do not influence 
the gravitational potential well. This gives an overview of the evolution of the system 
without costly computation. 





Figure 4: Parabolic passage of one massive companion (red) near another galaxy (blue) 
surrounded by 120 test particles (empty circles). The mass of the companion is one 
quarter of that of the primary. The numbers indicate the time in numerical units, with 
t = for the pericenter. (Adaptation of Figure 4 of Toomre & Toomre 1972.) 
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Figure 5: On the Earth, the gravitational potential of (mainly) the Moon deforms the 
shape of the oceans. According to the position of the Earth-Moon axis, we can experience 
high tide (left) or low tide (right), like in the bay of Fundy (between New Brunswick and 
Nova Scotia; photos by S. Wantman). 

Using several different configurations, Toomre & Toomre (1972) characterized the 
role of many parameters of the merger (like e.g. the mass ratio, the orientation of the 
spin of the progenitors...), and clearly showed that the perturbation of the gravitational 
potential of a disk by a point-mass intruder creates a bridge of test particles between the 
two galaxies and a tail in the opposite side of the disk, as shown in Figure 4. 

These features are very common in merging galaxies, and not only there... 



2 Tides 

2.1 On Earth as in Heaven 

Twice a day, a roll of water comes to the coast and makes the level of the sea go up 
and down (see Figure 5). Everyone knows that this physical phenomenon, called tide, 
is mainly due to the Moon but only a few understands it. The oceanic tides should 
be considered as a combination of many effects, including the attraction of the Moon, 
but also, the one from the Sun, the oceanic currents, the geometry of the coast and the 
bottom of the sea. Among them, the gravitational attraction of the Moon and the Sun 
plays a similar role to what occurs in galaxy mergers. 
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Figure 6: The attraction of the Moon (red arrows) is stronger in A than in E and C. The 
differential force (green arrows) is therefore opposite in A and in C. 

So, what happens? Let's concentrate on the Earth and the Moon. The gravitation (as 
described by Newton) between the Moon and an element of the seas acts proportionally 
to their masses (M and m) and the inverse square of their distance r: 



Consider five points (A to E) on the Earth (Figure 6). The force on a given mass element 
is stronger in A than at any other point while the one in C is the weakest. In the reference 
frame of the Earth, one has to subtract the force at the center of the Earth (in E). The 
differential force one obtains acts "outward" along the Earth-Moon axis (A and C), and 
"inward" along the perpendicular axis (B and D). That is why the seas are "pulled away" 
from the Earth in A and C, which makes for "high tide", while B and D are places of 
"low tides". 

This effect also applies to the atmosphere and the continents (Lindzen & Chapman 
1969). In the last case, the inertia and the friction make it almost invisible. Because the 
disturbance is less important than the binding energy of the Earth, the planet remains 
stable and binds its oceans, continents and atmosphere. Obviously, the Earth has a 
similar effect on the Moon. 

The same effect exists for much larger systems. In the case of galaxies, the material 
(stars, gas, dust ...) is bound by the gravitational potential of the galaxy itself. By 
modifying this potential, an external mass can disturb the equilibrium. The binding 
energy of galaxies is much smaller than that of planets, thus a perturbation like the 
close passage of an intruder galaxy can strip off their material over large distances. This 
explains the long plumes we observe in mergers. 



7 



General Introduction 



2.2 Prograde or retrograde, that is the question 

The mass ratio and pericenter distance of the two galaxies define the strength of the 
potential of the intruder with respect to the potential of the main galaxy. However, 
the tidal effect would be more efficient if it lasts for a long time. The duration of the 
perturbation corresponds to the time a mass element of the galaxy A sees the potential 
of the galaxy B, and thus is linked to the relative velocity of the two progenitors. 

For a given value of the velocity (which depends on the orbits of the galaxies), two 
cases can happen. In the reference frame of the galaxy A, the internal angular momentum 
(i.e. the spin) of B can be aligned or anti-aligned with its orbital momentum. In the 
first case, called prograde encounter, the relative velocity of the material of B close to 
the galaxy A is low, and thus, the tidal effect of A on B acts longer. A contrario, during 
a retrograde encounter, this velocity is higher and the effect shorter. 



■ Box 1 



In the reference frame of the galaxy A, the orbital motion of the galaxy B, 
which has a rotational velocity Q, is characterized by a linear velocity 



v rb = ?"ab e r x Q e z . 



(2) 



For simplicity, let's assume that the internal and orbital motions of B are 
coplanar, i.e. both rotational vectors are (anti-) aligned with e z . 



to 




.A 



Therefore, the vector of the internal rotation can either be +u>e z , or — ooe z . 
Thus, the velocity of an element of mass of the galaxy B, situated at a distance 
re relative to the center of its host galaxy is 



Vmt = (~tb e r ) x (±u e z ) 



(3) 
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and the norm of total relative velocity reads 

|v rb + v int | = r AB ^Tr B w, (4) 

where the sign of the second term depends on the alignment of uj with e z , i.e. 
with f2. For a prograde encounter, the spin and the orbital motion (fl) are 
coupled (i.e. aligned). Therefore, the relative velocity is lower (tab O-tbw) 
than for a retrograde encounter (tab + re oj). 



As the result of the encounter, the tidal structures would be much more pronounced 
in a prograde illustrated in Figure 7. 

Simple arguments allows to predict the global behavior of a pair of interacting galax- 
ies. However, the precise, fine description of the tidal field requires some mathematical 
definitions. 



2.3 Tidal tensor 

In the previous paragraphs, we have seen that tides arise from the gravitational at- 
traction of a body on a spatially extended object. Hence, it is natural to introduce the 
spatial derivative of the force to characterize them. Let's consider, for example, a small 
galaxy D ("dwarf") orbiting around a bigger, spherical one a noted MW. According to 
Newton's second theorem, we can consider MW as a point mass, situated at its center 
of mass. The net acceleration a of a star a of the dwarf is the sum of the external effect 
of MW at its position r a , and the one of all other stars in D: 

a a = a ext (r a ) + ^ a /3( r «)- (5) 

Because the center of mass of D undergoes only the acceleration of MW, we can reproduce 
the differential scheme we used for the Earth (Figure 6) by switching to the reference 
frame of the dwarf galaxy: 

a Q , D = [a ox t(r a ) - a cxt (r D )] + ^ a p{ r a)- (6) 

P + CL 

We note that the distance 6 = r a — td between a and the center of D is much smaller 
than the distance to MW. Therefore, one can linearize the previous expression and get 

a a , D = 8 <9a ext + ^ a^(r a ). (7) 

a Cosmological simulations show that this situation is extremely frequent and may be a fundamen- 
tal step of the creation of the galaxies, within the hierarchical scenario. In the nearby Universe, the 
Magellanic clouds orbiting around the Milky Way are an example. 
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Figure 7: Toomre-like restricted simulations of prograde-prograde (a), retrograde- 
prograde (b), prograde- retrograde (c) and retrograde- retrograde (d) encounters between 
a major galaxy (blue) and an intruder (red). The mass ratio is 1:4; each galaxy is made 
of 500 particles; the intruder is set on an eccentric orbit (black dotted line). Arrows 
indicate the spin of the two disks. The orbital angular momentum is the same in all con- 
figurations. Large tidal features originating from a disk are only visible in the prograde 
case while a retrograde encounter leaves the galaxy fairly unaffected. 
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Using the Einstein's summation convention, we have 

<D = 5^'(4 xt )+^4(r a ), (8) 

where the effect of the external potential on the dwarf galaxy is described by the term 

T ji = & « xt ) (9) 

which is the j, i term of the 3x3 matrix T, called tidal tensor. Knowing its components, 
it is possible to compute the effect of the tidal field at any position within the dwarf. 
Because a l ext can be derived from the gravitational potential — <9*0 ext , we note that the 
tensor is symmetric 

T ji = -^«9> ext = -<9^0 cxt = TK (10) 
Examples of the tidal tensors of classical density profiles are given in Appendix B. 



2.4 Introduction to compressive tides 

Aside from being symmetric, the tensor is real-valued and thus can be set in an 
orthogonal form. Three eigenvalues {Aj} measure the strength of the tidal field along the 
associated eigenvectors. The trace of the tensor reads 

Tr(T) = A * = = -V 2 ext . (11) 

i 

Poisson's law gives a link with the local density p (see Section A. 1.1, page 195): 

Tr(T) = -AirGp < 0. (12) 

This means that the trace of the tidal tensor cannot be positive, and thus that it is im- 
possible to compute simultaneously three positive eigenvalues. Remains the possibilities 
of two, one or no positive eigenvalues, as noted by Dekel, Devor, & Hetzroni (2003). In 
the Earth-Moon case (recall Figure 6), we saw that the differential forces act "outward" 
along the planet-satellite axis, while they are orientated "inward" in the other two di- 
rections, giving two negative A's. In analogy with solid bodies, the tides exert a stress 
along the Earth- Moon axis, and a compression along the remaining axes. As for those of 
the stress tensor, the eigenvalues of the tidal tensor denote a compression when negative 
and an extension when positive, along the associated eigenvector. 

Figure 8 sums up all the configurations, labelled by the number of negative eigenval- 
ues. Remember that the first case (0) is not possible in the Newtonian gravity, as it is 
linked to a negative density. Cases (1) and (2) are well known situations in, e.g. clus- 
ters of galaxies (Valluri 1993) or spiral waves (Gieles, Athanassoula, & Portegies Zwart 
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0< k l <K 2 <K 3 X 1 <0<X 2 <^3 A- 2 < <X 3 X 1 <K 2 <k 3 <0 



Figure 8: Different cases of tidal effects on an element (here, a green cube). In the con- 
figuration (0), all eigenvalues A« are positive and the tides act as an extensive effect (blue 
arrows). In the case (1) (respectively (2)), the tides are compressive along one (respec- 
tively two) direction(s) (red arrows). In the last case, the tides are fully compressive: all 
the eigenvalues are negative. 

2007). The last case, where all three eigenvalues are negative has been overlooked in 
the literature. In this configuration, the tides are fully compressive. In the rest of the 
present manuscript, we will refer to such a case as compressive tidal mode. We will come 
back to this important point in the next Chapters. 

Compressive or not, the mathematical representation of the tides allows us to describe 
the large scale effects on smaller scale objects, like dwarf galaxies or star clusters. 



3 Star clusters 

3.1 Collection of stars 

Any astronomer will tell you that clusters of galaxies contain galaxies, which contain 
clusters of stars. Indeed, as the galaxies are often organized in clusters, the stars that 
form the galaxies are themselves grouped in clusters. Discovered first during the 17th 
century, many star clusters have been catalogued by Charles Messier in 1774, without 
knowing that these objects were groups of stars. This update came later, in 1791 when 
William Herschel looked at M5 with his 40 feet-long telescope. He saw lots of stars in 
the cluster and even a core so dense that he could not resolve all its components. 

With more and more powerful observational techniques, it has been possible to detect 
numerous clusters and to resolve them into stars in the Milky Way and nearby galaxies 
(e.g. Andromeda). Two types of clusters revealed themselves (see Figure 9): 

• the globular clusters present a smooth spherical shape and a very high density in 
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Figure 9: Two star clusters in the constellation of Gemini. NGC 2158 (top left) is 
a globular cluster of ~ 1.5 Gyr while M35 (bottom right) is an open cluster 10 times 
younger. The blue massive stars in NGC 2158 have already exploded as supernovae 
when many of them are visible in M35. (Image by J.-C. Cuillandre, CFHT.) 

their center. They are often old (~ 1 — 10 Gyr) and very massive (~ 10 4 — 10 6 M Q ; 
see Meylan & Heggie 1997 for a review); 

• the open clusters are much younger than the globulars, more metal-rich, and have 
low masses, as for the Pleiades or the Hyades clusters (< 5000 M Q ). 

Despite their differences, the two classes of star clusters share similar sizes, which leads 
to a large spread in their central density (10 1 — 10 6 M pc -3 ). This historical classification 
mainly relies on observations of clusters of the Milky Way, a galaxy that does not undergo 
a strong burst of star formation. However, in galaxies showing high star formation rates, | SFR 
recent observations revealed the existence of young globular clusters, having the density 
and mass of globulars but a metallicity and mass function similar to those of open clusters 
(see e.g. Holtzman et al. 1992; Whitmore et al. 1993; O'Connell, Gallagher, & Hunter 
1994; Watson et al. 1996; Maiz-Apellaniz 2001). The high mass-to-light ratio for these 
young stellar populations makes the clusters detectable over large distances, that is, in 
other galaxies than our Milky Way. They are suspected to be at the formation stage of 
the classical globular clusters, as those we observe in our Galaxy, provided they survive to 
the possible destructive effects. Indeed, we know that most of the stars form in clusters 
(Lada & Lada 2003) but, in the Milky Way, only 1% of them are still members of a 
cluster. This suggests a strong destruction factor, from internal and external causes. We 
will further explore this point in this manuscript. 
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3.2 So many radii 

In the twentieth century, the growing body of observational data allowed a precise 
modeling of star clusters. The mass profile has been described by e.g. Plummer (1911), 
Jeans (1915), Michie (1963), King (1962) or King (1966), followed by evolutionary studies 
(see e.g. Spitzer 1940; Henon 1961). These theories have led to a large number of physical 
quantities, among which the radii, which may be confusing, even today. For the sake 
of clarity, the definitions of some of these quantities are given below, considering the 
example of a star cluster orbiting in a galaxy. 

• The tidal radius refers to the radius where the density of the cluster drops to zero. 
It is well defined for analytical profiles like the empirical King (1962) or theoretical 
King (1966) profiles, but very difficult to determine observationally. It is important 
to note that even the isolated objects have a tidal radius, which has stricto sensu 
no relation with tides. 

• The Jacobi limit corresponds to the boundary of the region inside which the cluster 
is gravitationally dominant. For simplicity, it is often defined as the distance from 
the center of the cluster to the first Lagrangian point. In this case, it might be 
called Jacobi radius (see e.g. Baumgardt et al. 2010), but also Roche sphere or Hill 
sphere. Rigorously however, the Jacobi limit is not spherical. Outside this limit, 
the tidal forces due to the galaxy strip the material off, thus creating a cut-off limit. 
Therefore, the Jacobi limit is necessarily larger than or equal to the tidal radius. 

Slightly different definitions also exist: the Roche limit usually applies to planets or 
satellites. It describes the limit where the tidal forces due to a larger object overcome 
the cohesion forces of the planet, thus destroying it. The Roche lobe refers to stars 
or binary stars where the tidal field is negligible but replaced by the thermodynamic 
pressure. When a star fills its Roche lobe, its material escapes from the gravitational 
bound. 

In the following, we will focus on the properties of a star cluster and thus, using the 
historical definitions, we will refer to the radius where the density drops to zero as the 
tidal radius, and to the radius depending on the tidal field of the host galaxy as the 
Jacobi limit. 

3.3 Formation and environment 

If the general process of the formation of a cluster is well known, many details remain 
unclear and thus, out of reach of numerical simulations. We do know that a cluster forms 
from a giant molecular cloud (> 10 4 M ) that cools down (to a temperature of < 10 K), 
fragments and then collapses. Therefore, the temperature, chemical composition, density 
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Figure 10: First steps in the life of a star cluster. Left: the first 130 stars of NGC 7129 
formed in a gaseous and dusty nursery which they begin to spread away, in a bubble 
shape. Middle: in the Small Magellanic Cloud, the young stars of NGC 346 have al- 
ready expelled a fraction of the surrounding gas but are still embedded in it. Right: in 
NGC 3603, the cluster has now swept out most of its gaseous content. (Images by T. 
Megeath, Spitzer; A. Nota, HST; J. Maiz Apellaniz, HST.) 



and pressure must control this process. In addition, shocks (Barnes 2004), turbulence 
(Elmegreen & Efremov 1997; Mac Low & Klessen 2004), ram pressure (Kapferer et al. 
2009) and magnetic fields (Verschueren 1990) are also regulators of the formation and 
must be correctly parametrized in the simulations. With so many knobs requiring a 
fine tuning, most of the studies (observational or theoretical, recent or older) focused 
on the formation of clusters from a global point of view, characterizing their location 
(Mihos, Bothun, & Richstone 1993), timescale (Elmegreen 2000), rate (Schmidt 1959; 
Kennicutt 1998), or efficiency (Rownd & Young 1999; di Matteo et al. 2007), in isolated 
or interacting galaxies. 

After these phases leading to a proto-cluster, the first massive stars form in a clump 
of gas which they rapidly (~ 1 Myr) blow away thanks to radiation and stellar winds. 
Note that all the gas from the molecular cloud is not converted into stars: the fraction 
of mass actually processed is called the star formation efficiency e: | SFE 

(13) 



Mi 



total 



(generally of the order of a few percent), which leaves a fraction (1 — e) in the gaseous 
phase in the young cluster. 

At this point, stars themselves are not directly visible, but via the re-emission of 
their ultraviolet light by the dust, in the infrared wavelengths. The process continues UV 
with the destruction of the dust by intense UV radiation (Draine & Salpeter 1979), the 
formation of lower mass stars and the explosions of the first supernovae (Castor, McCray, SN 
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& Weaver 1975). As a result, all the gas left-over from the formation of the stars is widely 
expelled after ~ 10 7 yr. If the cluster survives these steps and remains bound by its own 
gravitation despite this mass-loss (see Goodwin 1997; Geyer & Burkert 2001; Boily & 
Kroupa 2003a, 2003b; Baumgardt & Kroupa 2007), it can begin its life in a gas-free 
region. Figure 10 sums up these major phases. 



3.4 Observations and completeness 

The link between observational and theoretical studies is however biased when the 
detection limits of the fainter clusters are reached. Once a star cluster is formed, the 
stars will evolve so that the total luminosity decreases, as shown by McCrady, Gilbert, 
& Graham (2003), Weidner, Kroupa, & Larsen (2004) and summarized in Figure 11. 

This underlines two points that one must not forget when interpreting the observa- 
tions: 



Figure 11: Evolution of the visual magnitude of clusters of mass 10 6 M (blue) and 
5 x 10 5 M (red). Here, the 5 x 10 5 M cluster forms 15 Myr after the more massive 
one. (Data points from Weidner, Kroupa, & Larsen 2004.) 
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• The brightest cluster is not always the more massive one. Indeed, the 5 x 10 5 M 
cluster shown in Figure 11 yields a higher luminosity than the other one (twice 
more massive) for ~ 10 Myr. A constant mass-to-light ratio can only be assumed 
for clusters of similar ages. 

• Old clusters (for all masses), yield low luminosities and may not be detected by 
the instruments. When searching for statistical results on a population of clusters, 
one has to take the completeness limit of the sample into account and correct for 
the faint clusters. 



In the end, the determination of the mass and age functions of clusters within a galaxy 
may be extremely biased by these two points because the mass estimate can yield signif- 
icant errors, and the age distribution may miss the high ages (i.e. low My) end. We will 
see in Chapter II that extra caution must be taken when interpreting these observational 
estimates. 



Box 2 

The lifetime of a star cluster depends on the disrupting effects it may 
experience, both internal and external. In this way, the rate of death would 
depend on the initial mass of the cluster or in other words, its binding energy 
(Spitzer 1958; Baumgardt & Makino 2003). Observational data however, seem 
to be more controversial. Despite the difficulty to probe the low mass end of 
the mass function (because of the limit of detection), two empirical models 
have been proposed: 

• a mass dependent disruption model (MDD, Boutloukos & Lamers 2003) 
gives the disruption timescale tdis oc M 7 with 7 ~ 0.62. 

• a mass independent disruption model (MID, Fall, Chandar, & Whitmore 
2005). 

Both theoretical and numerical (iV-body) studies favor the MDD model (see 
e.g. Baumgardt & Makino 2003; Gieles & Baumgardt 2008). For instance, 
an external (extensive) tidal field truncates the cluster to the Jacobi radius, 
which can be linked to the half-mass radius (Ostriker, Spitzer, & Chevalier 
1972) and finally, to the mass (Larsen 2004). Therefore, the disruption of the 
cluster due to the tidal field is linked to the mass. 

The MID model has been introduced to explain the evolution of the mass 
function of the Antennae galaxies (that we will see in more details in the 
Section II. 1.3, page 58). Then, it has been extended to other galaxies like the 
Small Magellanic Cloud (Chandar, Fall, & Whitmore 2006) or M33 (Sarajedini 
& Mancone 2007). 
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However, the completeness problem in these galaxies biases the age distri- 
bution of the clusters toward a oc r _1 law (Gieles, Lamers, & Portegies Zwart 
2007), which forbids to conclude on the role of the mass in the disruption. 

Again, one must be careful when addressing statistical issues like the mass, 
luminosity or age functions, because the sample considered might not be com- 
plete. In any case, one should not forget that the Antennae widely differ from 
the other galaxies used in the studies of both models in the sense that, as a ma- 
jor merger, the environment of the clusters is very different from those existing 
in quiescent galaxies. This last point is explored further in the Chapter II. 

Note that both models assume a power law for the cluster initial mass 
function (CIMF): tf(Af) = dN/dM oc M^, with p ~ -2.0. However, a 
different CIMF has recently been proposed by Larsen (2009) and Gieles (2009), 
with the form of a Schechter (1976) function: 



where (3 ~ —2.0 is the usual slope of a power law CIMF and M c characterizes 
the exponential drop of the mass function (M c ~ 2 x 10 5 M for quiescent 
spiral galaxies, see Larsen 2009). Gieles (2009) discussed in detail the param- 
eters and the implications of such a function on the properties of the cluster 
populations. 



Nowadays, the knowledge of star clusters, and especially their early life, in our Galaxy 
or others (Harris 1991), appears to be the key to the understanding of as various phenom- 
ena as the formation of stars and the evolution of galaxies (see Brodie & Strader 2006 
for a recent review). Several authors explored the dynamics or the stellar populations of 
clusters to derive cause-to-effect relations between the clusters themselves, internal (e.g. 
stellar evolution in Hurley, Pols, & Tout 2000, mass segregation in White 1977 or mass 
loss in Hills 1980) and external (e.g. passage through the galactic disk or a spiral arm 
as in Gieles, Athanassoula, & Portegies Zwart 2007 or the perturbation due to the tidal 
field described in Lamers et al. 2005) processes. Recent literature suggests to link the 
formation of the clusters to their environment, i.e. major events in the life of their host 
galaxy. Indeed, clusters and their stellar populations can trace back the history of their 
host (see the inversion methods of Ocvirk et al. 2006), what is, in particular for mergers, 
a fundamental aspect of our understanding. Both galactic and cluster scales are coupled 
as one strongly influences the other through stellar winds, bubbles, chimneys, chemical 
enrichment or tides. The present document focuses on the last point. 




(14) 
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4 This thesis 

We have shown above that large-scale interactions between galaxies lead to the devel- 
opment of tides that act at cluster-scale. Putting all these points together, one naturally 
ends with the title of the present thesis: 

Dynamics of the Tidal Fields & Formation of Star Clusters in Galaxy Mergers. 

The work presented in this document tries to bring new elements to the giant jigsaw 
puzzle of interacting galaxies and stars clusters. Using numerical experiments, it aims 
to show that the large scale structures of galaxies give rise to phenomena that directly 
act on the smaller scales structures, like the star clusters. Among these effects, we 
concentrate on the tidal field in interacting galaxies and detail its bimodal nature by 
introducing and extensively developing the concept of compressive tides. This thesis 
suggests that the compressive tidal mode participates in the formation and early life 
of star clusters by preventing their dissolution and even helping their creation. This 
hypothesis is confirmed by an intensive comparison between our numerical results and 
the observations of a well-known pair of galaxies (the Antennae), and then extended to 
other systems. This document is organized as follow: 

• The first Chapter presents the numerical methods used during this work, to model 
galaxy mergers and to compute the effects of the tides. It first briefly explains 
the concept of iV-body simulations and details the algorithms used. The numerical 
technics and codes are then presented, as well as the test suite developed to validate 
the methods. 

• Their applications are shown in the second Chapter. There, we present the Anten- 
nae galaxies, from an observational point of view and then, the numerical model 
we created. We develop the results obtained about the tidal field of this model. 

• The strong link found between the existence of a compressive tidal mode and the 
formation of star clusters and tidal dwarf galaxies is detailed in the third Chapter, 
thanks to a comparisons with the physical quantities that monitor star formation. 

• Chapter IV extends these conclusions to a broader range of parameters, showing 
how robust the results are. Many cases of interacting galaxies are covered, with no 
particular link to observed systems. The Chapter particularly presents the influence 
of the distance between the galaxies during the interaction. 

• In Chapter V, we take one step backward and explores the link between the mass 
profile of the dark matter halos and the tidal field. We show that the choice of the 
profile strongly influences the character of the tidal field which in the end may, or 
not, support the formation of substructures. 
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• The numerous questions raised during this PhD and the perspectives they bring 
are presented in the Conclusions. 

• Finally, the Appendices bring more details to some points of the main text. The 
interested reader will find there an analytical complement (with many detailed 
derivations) that may help her/him to better understand this manuscript. 



Gas 
Tides 
Stars 
Gravity 
Clusters 
Galaxies 
Computers 

Now that all the characters have been introduced, 
let's begin the/to play. 
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Chapter 1 

Numerical techniques 



Overview 



This Chapter presents the numerical methods used in 
this work. It explains the choices made and the tests 
performed to validate the study. 



J 



1 TV-body simulations 

1.1 Smoothing the singularity 

In astrophysics, the iV-body gravitational problem usually consists in solving the 
equations of motion of a system of N particles using the Newtonian description of the 
gravitation. The acceleration of the i-th particle is simply given by 



where ext represents an optional external potential (which we will consider null here- 
after). Once the initial conditions in position and velocity are set, the solution can be 
found by solving the N systems of non-linear second order ordinary differential equations 




(1.1) 



< 




(1.2) 



dt ' 
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via e.g. an Euler or a leap-frog integration scheme. In the last case, the solutions in 
position %i and velocity Vi are updated with 



However Equation 1.1 yields a singularity when the distance \vj — r$| between two 
particles approaches zero. If such a situation occurs, an integration with a constant 
timestep dt would make the (a\ dt) term extremely large, leading to (unrealistic) high 
velocities and finally unbinding close particles. The obvious solution would be to use an 
adaptive timestep, that gets smaller when the distance between two particles becomes 
too short. All the other particles would then evolve with the timestep of the encounter 
(i.e. several order of magnitude smaller than usual). Because this situation may occur 
repeatedly in the system, this can considerably slow down the calculation. 

A more efficient solution, invented by Aarseth (1963), consists in modifying the law of 
gravitation on the very small scales, by introducing a smoothing (or softening) length 3 
e, in Equation 1.1 (with ext = 0) 



With this method, the singularity does not exist as the denominator never vanishes. 
This avoids using very small timesteps, but at the price of introducing an error in the 
force. Note that e can be used in other functional forms, called kernels, that change 
the expression of the approximated Newtonian force and can ensure, e.g., continuous 
derivatives (see Dehnen 2001 for a nice introduction, and the Section B.2, page 211 of 
the present document for some analytical expressions). 

Small values of e render a noisy system because of the finite number of bodies used. 
At the opposite, large values increase the bias in the expression of the force. The value 
of the smoothing parameter had traditionally been chosen in an ad hoc way for a long 
time, until that Merritt & Fridman (1996) and Athanassoula et al. (1998) proposed 
an objective method minimizing both sources of error. The optimal smoothing length 
decreases with the number N of particles, roughly as N~ 0Ai (fine details depending of 
the actual distribution of particles within the system). 

1.2 TV-body versus physical units 

In Figure 1.1, we show a typical output from an iV-body simulation. Can you tell if the 
system modeled is a globular cluster, an elliptical galaxy or even a cluster of galaxies? As 

a This approach is similar to the SPH method that smoothes the physical properties over a scalelength. 




(1.3) 




(1.4) 
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Figure 1.1: iV-body realization of a Plummer sphere (N = 5 000). The axes of this Figure 
do not have explicit units. Instead of physical units (like kpc), we have kept the iV-body 
units system. 



one rapidly notices, the axes of the figure do not have units like kpc, therefore answering 
the question is not possible. 

In the 1980 's, with a growing number of simulations, the iV-body unit system has 
been introduced (Heggie & Mathieu 1986) in order to allow for an easy comparison of 
the models from various authors. It sets 

G = M = -AE = 1, (1.5) 

where G is the gravitational constant, M is the total mass and E the energy. The -4 
factor for E corresponds to a virial radius set to unity (see Section E.2, page 235 for 
details). One can deduce the units of distance d and time t: 

GM 2 /T , 

d=~^ (1.6) 

d 3 



GM 



G VJi^- (L7) 

We note that this system of units is not physically constrained, as soon as a dimension- 
less quantity like GMt 2 /d 3 is left unchanged. Indeed, one can set the size of Figure 1.1 to 
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20 pc and its mass to 10 5 M Q . In that case, one of the 5 000 particles weights 20 M and 
one spatial unit represents 0.7 pc. This leads to a time unit of ~ 3 x 10 4 yr. However, by 
setting the values to bigger scales, e.g. a size of 30 kpc and a mass of 10 11 M Q , one gets 
a timescale of ~ 1.5 Myr. In the first case, Figure 1.1 would represent a globular cluster, 
while in the second, it would more likely be an elliptical galaxy. That is why purely 
gravitational simulations 5 are freely scalable, and thus must be correctly documented to 
avoid misinterpretations. 

In the modeling of galaxies, the stellar encounters are not relevant regarding to the 
lifetime of the object, and thus, such systems can be considered collisionless. 



Box 1.1 

Depending on the scales considered in a system, the role of each physical 
phenomenon may not be of the same importance. For instance, the dynamical 
timescale is linked to the typical size and velocity: 

R . . 

^dyn = — , (1-8) 
V 

and is approximately equal to the free-fall time: 



^"WaS? (L9) 

This must be compared to the relaxation timescale and to the lifetime of the 
object 



N 
"ln(iVy 



trelax ~ 0.1 t dyn (1. 10) 



(see Binney & Tremaine 1987). The following table gives order of magnitudes 
for the timescales of several systems: 



system 


*dyn [yr] 


^relax [y r ] 


globular cluster 


10 5 


10 9 


galaxy 


10 8 


10 16 


galaxy cluster 


10 9 


10 11 



In the case of a globular cluster, the relaxation time is comparable to the 
lifetime of the cluster and thus, stellar encounters should not be neglected. The 
(violent) relaxation processes influence the internal properties of the cluster 
via, e.g. mass segregation (White 1977). 



The introduction of hydrodynamical aspects like a cooling function adds a constraint on the 



timescale. 
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For galaxies however, the relaxation time is much larger than the typical 
lifetime (~ 10 10 yr), which indicates that they are collisionless systems in 
general. In the case of clusters of galaxies, both the relaxation time and the 
lifetime are comparable and thus, the conclusion on the collisional nature of 
such systems must take into account the details of their structure. 



The present Chapter describes the numerical methods and the associated tests, which 
are not linked to any physical object. Therefore, the numerical units will be used here, 
before switching to the physical system in the next Chapters. 

1.3 Barnes, Hut and Dehnen meet Newton under a (apple) tree 

For an entire system made of N particles, the direct method requires to compute 
(N — 1) forces on each of the N particles, i.e. N (N — 1) ~ N 2 calculations. For physical 
systems like galaxies which require a large number of particles, the N 2 scaling becomes 
the major limitation to the high resolution studies. Nowadays, the largest simulations 
counts up to 10 10 particles (see e.g. the Millennium-II simulation, Boylan-Kolchin et al. 
2009) and obviously take advantage of other techniques. 

Among many methods, we quickly describe here the tree code, which has been in- 
tensively used during this thesis. Barnes & Hut (1986) introduced an inventive way of 
accelerating the computation considerably. Their idea is to group particles in branches 
(or cells) and to consider only the interaction between leaves and branches. 

First, the space is divided into virtual cells according to the following recursive scheme: 
if a cell contains more than one particle, it is divided into daughter subcells of exactly 
the half width (see Figure 1.2, left). A cell with only one particle is called a leaf, while 
cells that have been split into subcells are branches. A pseudo (i.e. virtual) particle is 
then assigned to each branch, placed at its center of mass and contains the total mass of 
its cell. All the non-virtual particles in the system feel the effect of the pseudo-particles 
(instead of all the other real particles) if they are situated far enough (see Figure 1.2, 
right). 

A tolerance parameter 6 (~ 1) is used to define the distance criterion. A cell (branch) 
of a size £, at a distance D of a given particle (leaf) is considered as a group if £/D < 9. 
If not, its subcells are examined with the same criterion, and this recursively until (i) the 
criterion is fulfilled or, (ii) the subcell is a leaf (containing only one particle). (Note that 
the case 6 = corresponds to a direct summation code, as the criterion on 9 is never 
fulfilled.) 

The total number of interactions is thus considerably reduced. This tree code brings 
the scaling down to N log (N) , which accelerates the computation (and allows higher 
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resolutions). 

Few years later, Dehnen (2000) improved the tree code by considering cell-cell instead 
of particle-cell interactions. When building the tree, each particle "remembers" in which 
cell it belongs. Then, all the particles in the same cell get the same effect from external 
cells. Figure 1.3 describes the three main algorithms: the direct summation, the tree 
code and Dehnen's code. 

Although the actual gain of this method depends on the initial conditions of the system 
and on the rules used to build the tree, it has been noted that the scaling was reduced 
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Figure 1.2: Left: the entire space (here shown in 2D, for simplicity) is split into virtual 
cells and subcells. A (sub)cell is divided into daughter (sub)subcells if it contains more 
than one particle. These one-particle cells are called leaves. Right: then, the acceleration 
on a given particle (red dot) is computed thanks to an angular criterion 9 (see a schematic 
definition at the bottom-left corner of the right panel, and the text for details). If a group 
of cells (called branch, green squares on the figure) is viewed under a small angle from 
the red particle, then it is far enough and its effect can be approximated. All the leaves 
it contains are replaced by a virtual particle (green star) whose mass is simply the sum 
of the masses of all leaves, and situated in their center of mass. If a branch is too close 
to the red particles, then its leaves are left as it and considered separately, as in a direct 
code. In other words, if a group of particles is far and small enough from an external 
point, then the gravitational effect of these particles is taken into account as a whole, 
instead of separately, which speeds up the computation. 
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Figure 1.3: The direct summation (a) considers the interaction between all particles, 
whatever their respective distance is. The classical tree code (b) gathers particles at 
long distance into a single cell, while close particles remain treated as in the direct code. 
Dehnen's version of the tree code (c) takes advantage of the symmetry of the force and 
computes the interaction between the cells themselves. 



to iV (see Figure 1.4). Recent tests revealed that Dehnen's code (called gyrf alcON c ) 
running on a single processor, competes with implementations of classical iV-body codes 
distributed on up to heigh equivalent nodes (Fortin 2006). Because technical details of 
these codes are not in the scope of this manuscript, we refer the reader to their introducing 
papers and the references therein (Barnes & Hut 1986; Dehnen 2000; Dehnen 2001). 



1.4 Finding Nemo 

All the codes presented above are publicly available through a toolbox (Barnes et al. 
1988; Teuben 1995) called Nemo d that contains hundreds of smallish and bigger codes 
used in classical iV-body modeling. All these tools can be used in a shell, giving the 
parameters at the command-line or within an Unix shell script, as shown in Figure 1.5. 
It is therefore possible (and rather easy) to create two galaxy models, merge them into 
one snapshot, give them initial positions and velocities and make them collide, via e.g. 
gyrf alcON. Visualization and diagnostic tools are also available. 

This study made intensive use of the Nemo package, which has helped to concentrate 
on the parameter space and the physics. However, the numerical parameters used by the 
codes have been chosen after many tests and checks which ensured that approximations 
and errors were under control (see Section 1.1.6). 



c GalaxY simulatoR using the Force ALgorithm with Complexity 0(N). 
d http://bima.astro.umd.eu/nemo/ [S 1 
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Figure 1.4: Performances of the tree code gyrf alcON (red) plotted as CPU time per body 
versus N, and compared to the classical tree code (blue). (Data points from Dehnen 
2000.) 



#!/bin/tcsh 

mkexpdisk out=mydisk.gal nbody=50000 

gyrfalcON in = mydisk.gal out=mydisk_gyr.gal tstop=50 
step=0.1 eps=0.01 theta=0.5 kernel=0 Nlev=6 

exit(O) 



Figure 1.5: Example of a Nemo script. This tcsh script chains the commands to setup the 
problem and solve it. It first creates a model of an exponential disk with 50 000 particles 
thanks to the command mkexpdisk and saves it in the binary file mydisk.gal. Then, it 
integrates the orbits of the particles with Dehnen's code gyrfalcON, over 50/0.1 = 500 
timesteps (50 iV-body time units), giving the values set by the user for the parame- 
ters (eps, theta, kernel, Nlev), and finally stores the result in the mydisk_gyr .gal file. 
Obviously, external programs can be inserted into the script to perform additional tasks. 
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1.5 Visualization 



As mentioned before, Nemo does contain visualization tools. However, for the purpose 
of this thesis, it appeared that they would not be flexible and powerful enough to allow 
a deep investigation. Therefore, we developed a brand new code whose only task is to 
render the iV-body systems modeled. Its name, Trackpart, came from its very first aim: 
it has been designed to plot the trajectory of selected particles and to go back and forth 
in time, tracking them (see Figure 1.6). 

Trackpart has always been developed in the scope that external modules or function 
could easily been added. That is why we chose IDL as programming language, a well- 
known tool in the astrophysics community. Starting from scratch, we created routines 
dedicated to the 3D viewing (xyz-rotation with the mouse, zoom, translation) and render- 
ing devices that efficiently compute on-the-go all the quantities needed (column density, 
line-of-sight velocity, isocontours...), depending on the 3D orientation. A graphical user 
interface has been built to read the input files, create images and movies quickly and 
easily. Most of the images of iV-body simulations in this manuscript have been prepared 




Figure 1.6: The Trackpart interface allows to 3D-rotate an A-body model, to zoom-in 
or -out on details, or select a region of interest, tag the particles in it and follow them 
along time (as done for the blue part of the tidal tail in this screenshot) or plot the orbits 
of particles (orange and green here). Other visualization modes (velocity field, column 
density ...) are also available. 
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using this tool. 

1.6 Testing the gravity 

In iV-body simulations, the most accurate result will come from a direct summation 
code with a very small timestep. As we have seen above, this would be an extremely long 
computation and is thus compensated thanks to approximations. Obviously, both these 
slight changes in the physical behavior of the system and the numerical treatment intro- 
duce errors in the final result. Hence, the parameters tuning the approximations must 
be chosen carefully, to get the best combination between accuracy and speed (keeping in 
mind that the solution will always be an arbitrary choice and thus depends on what is 
important for the user). 

First, we compare the results given by gyrf alcON with those of a classical tree code 
(called hackcodel). For this experiment, we place ourselves in a critical situation: the 
face-on collision of two disks. Let's consider two parallel identical exponential disks, each 
made of N = 50 000 particles, separated by a distance of five times their characteristic ra- 
dius. Without initial velocity, their mutual gravitational attraction makes them collapse 
on each other, as plotted on Figure 1.7. We expect the integrator to make the largest 
error on the computation of the accelerations when the disks go from quasi-isolated stage 
to violent collision, in other words, for a sharp evolution of the density. By quantifying 
the order of magnitude of the errors made in this case, one can evaluate the quality of 
the computation. 

We choose to monitor the thickness of one of the disks. For each snapshot, we extract 
the z-coordinate of every particle and compute the standard deviation a z . Results are 
shown on Figure 1.8. The long distance interaction first warps the disk (t < 9, panels a b 
c and d on Figure 1.7). At t ~ 10 and later, the thickening is accelerated by the collision 
which unbinds a fraction of the particles, resulting in a increase for a z . 

The two tree codes give comparable results (with relative differences of less than 
~ 5%). However, this does not mean that their results are physically correct. To check 
this, we monitor the conservation of energy and angular momentum while tuning the 
parameters of gyrf alcON (see Figure 1.9), like the smoothing length e, the softening 
kernel or the tolerance parameter 9. 

Except for a brief moment during the collision phases, the energy is well-conserved 6 . 
When tuning the parameters, we noted that the accuracy remains at a similar level for 
9 G [0.4,0.6], but the computation is sped up for large values of 9. This confirms that 
the default set of parameters (9 = 0.4, kernel PI, see Equation B.49, page 212) allows 
a correct usage of the code, at this resolution (i.e. N ~ 10 4 — 10 5 ). For higher number 

c Note that the angular momentum is de facto conserved by the intrinsic symmetry of the tree built 
by gyrf alcON. 
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Figure 1.7: Major steps in the evolution of a face-on collision between two exponential 
disks. Initially (panel a), the disks are separated by 5 times their characteristic radius 
and without initial velocity. As the gravity gathers them (panels b, c and d), the central 
region of a disk goes deeper into the potential well of the other, while the external parts 
are less rapidly affected. After the first collision (panel e), the disks are dynamically 
heated and become thicker. The remnants slow down considerably (panels f and g) and 
fall back on each other (panel h) before being indiscernible. The panels corresponds to 
the computational times 1.0, 8.0, 8.5, 9.0, 9.5, 10.0, 10.5, 11.0 and 18.0. 
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of particles however, e must be changed, according to the relation of Merritt & Fridman 
(1996, see also Section 1. 1.1 of the present document). 

The errors due to the collision are mainly due to the finite timestep used in the leap- 
frog scheme: when the acceleration becomes large, the potential energy is not properly 
converted into kinetic energy and the total energy E is not exactly conserved, as visible 
in the upper panel of Figure 1.9. In other words, the code is too slow in updating the 
physical quantities, with respect to the underlying physics. This effect can be reduced 
by using individual adaptive timesteps. In this case, when a particle satisfies a certain 
criterion (e.g. it lies in a region of high density), its proper timestep is reduced so that 
time sampling becomes a less critical aspect of the integration, gyrf alcON offers to play 
with adaptive timestep. After some rapid tests, it appeared that this was not critical 
and that the default timestep scheme (with small errors but rapid) was sufficient for our 
usage. 

The broad conclusions of these tests are that (i) the smoothing length e has to be cho- 



b 




5 10 15 20 

time 



Figure 1.8: Evolution of the thickness of one of the disks during the face-on collision 
shown in Figure 1.7. The thickness is measured as the standard deviation a z of the z- 
coordinates of the particles. Both a classical tree code (hackcodel) and Dehnen's code 
(gyrf alcON) show similar results, within an error range of ~ 5%. The increasing of a z 
is well explained by the collision (see text for details). 
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time 

Figure 1.9: Variations of the total energy E = K + Q (with K and Q, the kinetic and 
gravitational energies), angular momentum L and virial ratio v = —2K/Q during a face- 
on collision, as plotted on Figure 1.7. We note that the system is not virialized initially, 
but only after t ~ 2. The energy and the angular momentum are well conserved since 
the beginning of the simulation, showing that the code is well parametrized. However, 
the two collisions are clearly visible in the energy plot (t w 9.2 and t w 10.5, panels 
e and g on Figure 1.7). Here, the code faces the problem of the finite timestep, and is 
forced to sacrifice the accuracy to keep the computation fast. Anywhere else, the energy 
is conserved, with an relative error smaller than 10 -3 . We also note that the simulation 
shows the same level of energy and virial ratio before and after the collision. The virial 
ratio shows some fluctuations, but always less than 10%, which reflects the physical 
evolution of the system. 
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Figure 1. 10: Cube used in the computation of the tidal tensor. The central point P is one 
of the particles of the A-body simulation, while the six points A are pseudo-particles, 
added to compute the tensor. 

sen according to the number of particles, (ii) the tolerance parameter 9 can be set ~ 0.4, 
and (iii) the kernel and the adaptive timestep are second-order parameters regarding to 
our study. With its very high efficiency, the code gyrf alcON is perfectly suited for the 
simulations of interacting galaxies of this thesis. 



2 Computing the tidal tensor 

Once the iV-body simulation is done, we have a proper description of the mass distri- 
bution, i.e. the gravitational potential. The next step of our study consists in evaluating 
the tidal tensor at various points of the simulation, to see how the tidal field evolves in 
space and time. 



To compute the tidal tensor in our simulations, we start with the positions and masses 
of the particles at a given time, as input data. Let's consider a particle P from our 
simulation. According to Equation 9 (page 11), we simply need to compute the terms 
^( a ext)- To do so, a pseudo-particle A (i.e. mass-less) is placed on each side of a cube of 
size 25 centered on P, as shown in Figure 1. 10. Then, the net acceleration from all the 
iV particles in the simulation is computed at the positions of the pseudo-particles. 

We note that, in this case, we have not retrieved a* xt yet, but the total acceleration. 
Thus, we need to subtract the "internal" effects (i.e. the acceleration due to P at the 
positions of the A particles), that we do explicitly knowing the kernel and the distance 
S between P and A. Once done, we pair the pseudo-particles of two opposite sides of the 



2.1 Method 



38 



2 - Computing the tidal tensor 



cube and perform the first-order differencing. For example, for the x-axis 

rpxi _ Q ext(^+) ~ a Lct(Ar-) /j -. ,n 

After doing this along the three directions, we get the nine components of the tidal 
tensor, in the base {e^e^e^} which is arbitrary. Therefore, the next step is to diag- 
onalize the tensor and get the three eigenvalues {Aj} and the associated eigenvectors 

M- 

■ Box 1.2 



Technically, the code first reads a snapshot (called source) from the simu- 
lation to get the positions and masses of all P particles. Then, it builds a new 
snapshot, called sink, containing only the pseudo-particles (A). Both source 
and sink are given to getgravity, a Nemo tool that uses the same libraries as 
gyrf alcON, with the same parameters as gyrf alcON which has integrated the 
orbits in the simulation. Then, getgravity retrieves the three components 
of the acceleration for all sink particles. Knowing 5, e, the kernel used in the 
functional form of the gravity and the mass of the particle P, the effect of the 
central particle P on A is explicitly computed and subtracted from the total 
acceleration. The first-order difference scheme is then applied to retrieve only 
the upper part of the tensor (which is symmetrical). These six components 
are passed to the direct method of Kopp (2008) that solves the characteristic 
equation of the tensor and derives efficiently and quickly the eigenelements. 
Note that the errors introduced by this method are generally close to machine 
precision and thus, well below these due to the computation of the acceler- 
ations. The experiment is then repeated for all snapshots. The complete 
algorithm is schematically summarized on Figure 1.11. 

A slightly modified version of this algorithm has also been implemented: 
instead of centering the cubes on particles of the simulation, they are placed 
on a regular 3D grid, independent of the mass distribution of the iV-body 
model. In this case, there is no need to subtract the effect of a central particle 
in the cube. 



In the end, while the gravitational potential is resolved over the length e, the resolution 
of the tidal field is S. 



2.2 Tests 

The only parameter of this method, 5, must be chosen carefully to minimize errors. 
This problem is similar to the choice of the optimal smoothing length e. When setting 
the value of S, one has to cope with two sources of errors: 
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• Large values: the cube tends to include all the particles that build the external 
potential. The approximation of the distant source (see Equation 7, page 9) is not 
valid. 

• Small values: the iV-body rendition of the mass distribution fails to resolve the 
fine structure of the potential. In addition, round-off errors in the calculation 
of the finite-difference are introduced when the term [al xt (A x+ ) — a l ext (A x -)] in 
Equation 1. 11 becomes close to the numerical precision. 

Therefore, we immediately have a lower limit, set by the accuracy criterion on the 
potential, with 5 > e. Practically, the round-off errors only kick in when S is much lower 
than e, meaning that this will not influence our choice. But the problem of finding a 
suitable upper limit remains. To address this question, we ran tests on iV-body renditions 
of well-known potentials and compared the results obtained with different values of 5 to 
the analytical solutions derived in Section B.l (page 205). 



Tidal tensor computation 

A'-body simulation 



creation of a model 
T 

orbits integration 
(gyrfalcON) 



gravity parameters 
(e,6, kernel) 



t 



Figure 1. 11: Schematic algorithm of the computation of the eigenelements {Aj, z/j} from 
the creation of the iV-body model. Steps represented with dashed lines are skipped in 
the grid-based version of the code. 



extraction of a snapshot 



extraction of m, x, y, z 



creation of the sink 



gravity calculation 
(getgravity) 



I 
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subtraction of the effect of P 



first-order difference 



diagonalization 
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The simplest case is the well-known Plummer model whose potential is 

= ~ ^— 2 . (1.12) 

V r o + r 

Here, M is the total mass and ro the core radius. However, this profile yields an important 
issue for its iV-body rendition: it does not have a truncation radius. Indeed, its finite 
mass spreads over the entire space, which is unrealistic. That is why the models we 
use (since Figure 1.1) are artificially truncated. They are set up by the tool mkplummer 
(within the Nemo package) which considers that 99.9% of the total mass is enclosed within 
a radius of 22.8042468. Knowing from Equation B.5 (page 206) that 99.9% of the mass 
corresponds to a radius of ~ 38.7ro, we deduce that ro ~ 0.589 and can now normalize 
the results to compare them with the theory. 

Let's begin with the iV-body rendition of the potential. First, gyrf alcON computes 
the gravitational potential at the positions of the N = 500 000 particles of our Plummer 
model, which we compare to the theoretical value. Results are shown in Figure 1.12. 

Note that the standard deviation between the numerical and the theoretical values is 
cr^ ~ 10~ 3 . The relative error is everywhere lower that 0.6% and thus, the potential is 
very well rendered. Next is the acceleration: 

\<r)\ = GM ' 3/2 - (1-13) 
(Tq + r 2 ) ' 

Figure 1.13 plots the norm of the acceleration at the positions of the particles and, once 
again, the comparison with theoretical value. In this case, one has a a ~ 3 x 10~ 3 ~ 3cr^ 
and the relative errors reach 5% locally. With the N = 500 000 particles of our test 
case, we have set e = 0.01. Note that the passage from the potential to the acceleration 
increases the level of error by one order of magnitude. 

Although the same experiment is possible for the tidal tensor, the analytical derivation 
of its eigenvalues is much more involved f . Instead, we concentrate on the T xx term of 
the tensor, for y = z = 0, i.e. through the center of the model: 

T xx (r) = -GM^ ^ (1.14) 

(tq + r 2 ) 

(see Section B.l.l, page 205 for a detailed derivation). Obviously, one will not find 
particles with y = z = in the simulation, so we need to use the grid-based version of 
the code (recall Box 1.2) to compute the tensor along the x-axis. At a given resolution 
(i.e. for a certain number of particles), we vary the size of the cube 5, using the same 
parameters to compute the gravity as in gyrf alcON. Figure 1.14 and Figure 1.15 show 
the results. 



f The determinant of the tidal tensor, which is the product of its eigenvalues could also be used but 
again, its analytical derivation may be long, even for simple cases. 
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2 4 6 8 10 12 14 



r 

Figure 1.12: Top: Measured (blue) and theoretical (red) potential of an iV-body realiza- 
tion of a Plummer sphere (N = 500 000). Bottom: Relative error between the potential 
measured and the theoretical value. The red curve corresponds to the mean error per 
radial bin. For clarity, only 10 000 points are plotted in this figure. 



Table 1.1: o~t for several 8's 



5 
0.01 
0.02 
0.05 
0.10 



ax 



xlQ- 

3.3 

2.8 

2.2 

1.8 
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Figure 1.13: Same as Figure 1.12, but for the norm of the acceleration. 



The standard deviations are given in Table 1.1. As for the passage from the potential 
to the acceleration, the computation of the tidal tensor from the acceleration increases 
the level of error by one order of magnitude. This means that our method used to 
compute the tensor has a quality comparable to those of the tree code, used to get the 
acceleration. 

Figure 1.14 exhibits small, yet significant errors which are mainly visible for r « 2.1, 
2.5, 3.4 and 4.3. The origin of these errors is a problem which is discussed in details in 
Appendix C. 

As predicted, large cubes smooth the potential and thus the terms of the tidal tensor, 
while smaller cubes do not correctly reproduce the fine structures. However, with a too 
large cube, the finite-difference approximation is not valid anymore and an error is made 
in our calculation of the tensor. As visible in Figure 1.15, the curves corresponding to 
large 5 are smooth and regular but quite far from the analytical value. A contrario, the 
small values yield noisy curves, but on average closer to the theoretical solution^. 



g Note that, while the noise comes from the finite differences scheme, the systematic errors (i.e. the 
noiseless signal) are mainly due to the approximation made on the value of the force, when using a 
kernel. 
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Once again, the choice is arbitrary. An intermediate value like 5 = 0.05 seems to take 
into account both sources of error, reasonably. As mentioned before, the size of the cube 
strongly depends on the local description of the potential, that is the softening length e. 
In the previous example, we had e = 0.01 and 5 = 0.05. This suggests 6 = Be. 

In principle, the mean density in a cube must be independent of the resolution, i.e. 
iV<5 3 must be constant. With this assumption and the relation for e from Merritt & 
Fridman (1996), it is possible to test the validity of our empirical relation 5 = 5e, as 
shown in Figure 1.16. 

We note that the empirical relation matches those from the density of the cubes over 
a couple of decades, but the deviation becomes too large for small numbers of particles. 
The differences introduced when using one or the other relation are comparable to the 
changes for 5, already discussed in Figure 1.14 and Figure 1.15. 




Figure 1.14: Top: T xx term of the tidal tensor of a Plummer sphere. The tensor is 
computed at 1 000 positions along the radial axis. Bottom: Relative error between 
numerical and theoretical values. Note that the largest values correspond to extremely 
small absolute values of the tensor, and thus are not directly relevant as a measurement 
of quality. 
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One must not forget that all scaling laws above depend on the actual distribution 
of the mass. Therefore, one should not rely too strongly on them and introduce some 
flexibility Obviously the precise determination of the optimal values of e and 5 is out- 
of-reach for non analytical systems. That is why, knowing that the values chosen would 
rarely be optimal, we decide to simply use the quantities (iVe 3 ) and (N5 3 ) as constants, 
in the range of our typical resolutions (~ 10 5-6 particles). 



2.3 Alternative method 



In addition to the finite differences method presented above, another simple approach 
is possible. In an iV-body system, the gravitational potential at the point of coordinates 
(x, y, z) is given by the sum of individual contributions of all the particles: 



A' 



(z, y, z) = ^2 Vki.x, y, z), (1.15) 



k=l 




Figure 1.15: Zoom-in on the inner region of Figure 1.14. Large values of 5 are generally 
smoother (less noisy) than small values but further from the analytical solution (red 
curve). 
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where is the potential of the k-th particle, evaluated at (x, y, z). Thus, the components 
of the tidal tensor at this position read 

T i *(x,y,z) = -d i d i <l>(x,y,z) 

N 



k=l 

N 

k=l 
N 

y,z), (1.16) 



k=l 



where is the tidal tensor of the k-th particle. Here, T u can be interpreted as the sum 
of the individual contributions of the particles to the tidal tensor. 

A simple approach uses a point-mass potential for the particles, which leads to a tidal 




Figure 1.16: The relation giving the optimal value for the softening length e (blue curve) 
from Merritt & Fridman (1996) is compared to those assuming a constant density in the 
cubes (red), and to the empirical scaling law (red dashed). The proportionality factors 
have been chosen to get the same value for 5 at a resolution of iV = 500 000 particles. 
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tensor of the form 

T«(x,y,z) = ™(3x i x j -S*^). (1.17) 

However, as we have seen above, a kernel is often used to smooth the potential of the 
particles (see Section B.2, page 211). Once the tidal tensor associated with the right 
kernel is known, one has to sum the contribution of all the particles, at the position 
(x, y, z). Finally, the total tidal tensor T can be set in diagonal form, as before. 

The obvious advantage of this method is to eliminate the errors due to the finite dif- 
ferences scheme, as presented in the last Section. However, the previous implementation 
makes intensive use of Dehnen's tree code to compute the accelerations on the cubes. 
Therefore, it can be applied efficiently to high resolution simulations, for a small compu- 
tational cost. It is also possible to implement a tree method for our alternative approach. 
However, the limited time of this thesis did not allow such an exploration and we have 
only implemented a (time-wise costly) direct method. 

The results for our Plummer test case are shown in Figure 1.17. We also illustrate 
here the effect of the kernel. The finite difference scheme (which uses a kernel PI, see 
Equation B.49, page 212) is a noisy version (a = 0.048) of the associated alternative 
method (cr = 0.040), where maximum errors are found where the gradient of the accel- 
eration is important (see Figure 1.13). However, the computational cost of the direct, 
alternative method does not worth the gain in precision for the solution. 

The results presented in the following Chapters have all be obtained with the first- 
order differencing scheme. 
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Figure 1.17: Comparison between the finite differences methods and the alternative 
method, for again a Plummer sphere with N = 500 000. The black curve represents 
the tidal tensor computed with a Plummer kernel for the particles, while the blue curve 
denotes a PI kernel (see text). The green curve is the result found with finite differences 
(kernel PI and 5 = 0.05). 
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Overview 

This Chapter aims to introduce the galaxy merger 
which has been used as an observational reference 
throughout this work. Then, it details the results pre- 
sented in Renaud et al. (2008), and in the Sections 4 
and 6 of Renaud et al. (2009). 



1 A traditional merger 

1.1 A long time ago, in two galaxies (not so?) far far away... 

The main difference between astrophysics and the other fields of science is that any 
experiment is impossible to realize. Therefore, our theories rely on the simulations that 
we constantly compare with the observations of as many objects as possible. To study 
a phenomenon (galaxy mergers in our case), one has to collect a lot of high precision 
data to validate the theories. That is why the most observed objects are often the most 
reproduced numerically. 

For mergers, the natural candidate is the Antennae galaxies (NGC 4038/39, 
Figure II. 1), named after the peculiar shape of their long tails, already observed in pho- 
tographic plates by Duncan (1923) a who noted "two faint extensions, like antennae, [...] 
both concave toward the west" (see also Schweizer 1978). Indeed, the Antennae are the 
closest major merger observed, and thus one of the most studied 5 . 

The distance to this system has been subject to discussions over the last decades. 
a Back then, the two galaxies were still considered as a single nebula. 

b During the 30 months of this PhD, the ADS system has referenced over 170 publications mentioning 
the Antennae, i.e. one every 5 days! 
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Such an estimate is of paramount importance as it is used to derive almost all the other 
physical quantities (like the age, the size and even the mass). Many studies established 
distances ranging from 6 Mpc to 30 Mpc (see e.g. Rubin, Ford, & D'Odorico 1970; 
Whitmore & Schweizer 1995; Fabbiano, Zezas, & Murray 2001; Hibbard et al. 2001; 
Zezas & Fabbiano 2002; Saviane, Hibbard, & Rich 2004). This controversy is mainly due 
to the large errors introduced by every method. 

Up to recently, most of the authors agreed on a value of 19.2 Mpc, derived from the 
recession velocity (Whitmore et al. 1999). Saviane et al. (2008) used photometry of 
the tip of the red giant branch to determine a distance of 13.3 ± 1.0 Mpc, much shorter 
than the commonly used value. Few months before (on December 18, 2007), Drake et 
al. (2007) had detected a new type-la SN (called SN 2007sr) in the southern tail of the 
Antennae. Schweizer et al. (2008) used it as a standard candle to find a distance of 
22.3 ± 2.8 Mpc, in agreement with models of large scale flows of galaxies (Tonry et al. 
2000). Clearly, the debate is still open. 

In the following, we adopt the common value of 19.2 Mpc, which gives a total extension 




Figure II. 1: The Antennae show two extended tidal tails as well as a complex web of 
structures in the central region. North is up, east is left. (Image by R. Gendler.) 
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of the Antennae (from the tip of one arm to the one of the other) of ~ 110 kpc. 
■ Box II.l 



If one assumes a different distance for the Antennae, the intrinsic scales will 
change. Let D = d'/d be the ratio between a distance d' (e.g. 13.3 Mpc) and 
the value adopted d (here 19.2 Mpc). Starting from the observational data, 
the apparent diameter and the observed magnitude of an object must be the 
same for all distances. This implies that the other physical quantities (e.g. 
size, mass, timescale) are different. The invariance of the apparent diameter 
gives the change in length: 

R' = RD. (II.l) 

Usually, the absorption term is implicitly included in the apparent magnitude. 
Here, we explicitly write the absorption A* as an additional factor to the 
magnitude without absorption, m*, to get the observed apparent magnitude: 

m obs = + A*, (II. 2) 

i.e. the quantity measured by the telescope. With the definition of the distance 
modulus, one can link m obs to the absolute magnitude M* 

m obs = M* + A* + 5 log ( — j . (II.3) 

m Q bs is also invariant, therefore 

^-obs = ^obs 

M* - Mi = ( A[ - A*) + 5 log (D). (II A) 

The term A+ includes all the effects of absorption, i.e. from the atmosphere 
(for ground based telescopes), the Milky- Way, the intergalactic medium and 
the inner absorption of the object itself. When varying the distance to the 
Antennae, only the intergalactic contribution will change. For instance, when 
going from 13 Mpc to 19 Mpc, the column density of the intergalactic medium 
increases by 6 Mpc. However, the densities measured there are so low that 
one can reasonably simplify the problem by assuming that the difference is 
negligible (Schmidt 1975 derived the order of magnitude of 10 -3 mag Mpc -1 ), 
in other words A+ — A+ « 0. 

The absolute magnitude is also linked to the luminosity L± through 

= —2.5 log {Li) + k, (II.5) 

where k is a constant. 
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The mass-to-light ratio T = Mj does not vary with the distance d, thus 
we have 

M, - .1/' -2.5 log f ^ J + k + 2.5 log ( ^- J - fc 



'M'\ 

51og(D) = 2.5 log I — 



M' — M D 2 , (11.6) 

which is the equivalent relation for the mass. Then, for the density 
(p oc M/R 3 ) and the timescale (t oc p -1 / 2 ), we get 

p' = P D- 1 (11.7) 

t' = t D 1 ' 2 . (II.8) 



In other words, if d' < d, the density of a cluster is higher if the Antennae 
are situated at d! than if they are at d. 

It is usual (see Box II. 2, below) to consider that the stars are formed when 
the density of the interstellar medium overcomes a threshold p c , which must be 
independent of the distance between the galaxies and the observer. Therefore, 
p c is reached more easily in dense galaxies, or in our case, when the Antennae 
are close. That is, if the Antennae are situated at a distance d' = 13.3 Mpc 
instead of d = 19.2 Mpc, the formation of stars must be much more efficient. 

Indeed, consider that we want to form a cluster of mass M (the same for 
all distances), from a molecular cloud of size l\ (respectively £[) and density 
p (respectively p'). We need to condense the cloud from its original size l\ to 
a smaller object of size £2 and density 



Pc 


= c p 


M 


M 






£2 


= c- 1/3 e 



where c is a constant. The laws of physics describing the formation of stars 
are independent of the distance and therefore 

Pc = Pc 
I I 

c p = c p 

d = cD. (11.10) 
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The mass being also independent of d, we have 






P '£[ 3 = P i\ 
t\ = l x D 1/3 , 


(11.11) 


which finally leads to 








£' 2 = 4- 


(11.12) 


Once we have these relations, let's have a look at the influence of d on the 
strength of the tidal field denoted by A, which has the dimension of Gp, that 
is 




A' = A D-\ 


(11.13) 


The tidal energy per mass element for the cloud of size l\ is 






-^tidcs,! = — -Aa i\ 


(11.14) 


(see Section E.1.3, page 


234), which means that 






-^tidcs.l = #tides,l D 1/3 . 


(11.15) 


Once the stars have formed, the region to consider has the 
tidal energy becomes 


size £2 and the 




-^tidcs,2 — ^tides,2 D . 


(11.16) 


When compared to the 


gravitational energy Q per unit mass, 


we have 




-^tides,l -^tides,l 


(11.17) 


for the cloud, and 








-^tides,2 -Etides,2 


(11.18) 


for the cluster. 






In the end, the role of the tidal field remains the same or is slightly under- 
estimated if the distance to the Antennae drops from 19.2 Mpc to 13.3 Mpc. 
Note, however, that all order of magnitude computations are not affected by 
these uncertainties. 
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1.2 Meet a strange animal 

During a symposium in (formerly West-) Germany, Schweizer (1978) referred to the 
Antennae as an "animal from the zoo of Alar Toomre" . He categorized the galaxies in 
the subspecie of the "animals which always huddle in pairs, each individual having one 
long tail" . Even if the long tidal tails are the first visual signature of the Antennae, most 
of the high and very high resolution observations have been performed in the central part 
of these galaxies. They revealed an incredibly rich and complex web of substructures, a 
sign of a strong activity in the last ~ 10 7 yr (Figure II. 2). 

Velocity maps and early simulations of the Antennae suggest that the two nuclei have 
already undergone a passage next to each other and that we are observing the second 




Figure II. 2: The central part of the Antennae gathers the two nuclei of the progenitor 
galaxies as well as numerous star clusters created during their encounter. The area 
between the two nuclei filled with dust (brown) and molecular clouds (red/pink) is called 
the overlap region. Note that the tidal tails (not shown here) crosses, so that the southern 
tail is associated to the northern nucleus (NGC 4038) while the northern smaller arm is 
linked to the southern core (NGC 4039). (Image by B. Whitmore, HST.) 
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Right Ascension (J2000) 



Figure II. 3: Intensity-weighted Hi velocity map of the Antennae (left) and their central 
region (right). The color scale is the same for both images, with green indicating velocities 
close to systemic. On the right-hand side panel, one can distinguish the rotation of the 
nuclei, especially for NGC 4039, in the bottom part. (Adaptation of Figures 2 and 9 of 
Hibbard et al. 2001.) 



one c . These two passages would be responsible for the fireworks we see today. One 
might notice the major features of this region: the arc-shape structure in the northern 
part surrounding the nucleus of NGC 4038, and the overlap region between the two 
cores. In both lot of star clusters are visible. It is likely that they originate from 

the tidal interaction, because they show ages much younger than those of the galaxies. 
This recent and violent episode of star formation, called starburst, would be typical for 
interacting galaxies. However, the exact recipe of this formation is still subject to debate, 
because of the complexity of such a system. 

In order to better understand the Antennae, one has to gather morphological and 
kinematical data, from the numerous different observational sources. The work of 
Hibbard et al. (2001) brings many clues, through Hi maps. E.g., in Figure II. 3, one 
can easily distinguish the two progenitors thanks to steep velocity gradients. Note also 
the peculiar distribution of Hi in the northern tail, much shorter than the southern. This 
suggests that the major merger is not a symmetric object, and that NGC 4038 (northern 
nucleus, southern tail) dominates the present distribution of matter. 

c Note however that this contradicts the classification of the Antennae as the first step of the Toomre 
sequence (see Figure 2, page 3) while the next step (represented by the Mice) clearly shows an undergoing 
first passage. See also Karl et al. (2010). 
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The merger has been observed in all the other wavelengths: e.g. radio (Hummel & 
van der Hulst 1986), far-IR (Bushouse, Telesco, & Werner 1998), mid-IR (Wang et al. 
2004), near-IR (Brandl et al. 2005), UV (Hibbard et al. 2005) and X-ray (Baldi et al. 
2006). With such a mine of informations available, drawing the portrait of the Antennae 
is out of the scope of the present document. Important facts about the star clusters and 
TDG | the tidal dwarf galaxies are quickly presented in the following. However, the reader is 
strongly encouraged to consult these papers and the references therein for more details. 



1.3 Observations of star clusters... 

Over the last two decades, both observational and theoretical studies have shown 
that intense episodes of star formation take place in interacting galaxies. As the closest 
sample, the star clusters of the Antennae have been explored by several authors who have 
detected some 1000 of them in the central region of the merger (Whitmore & Schweizer 
1995; Meurer 1995; Mengel et al. 2005; Whitmore, Chandar, & Fall 2007; Bastian et al. 
2009; Whitmore et al. 2010). 

Using deep images in the V band, Zhang & Fall (1999) derived the mass function of 
these clusters and found an important discrepancy with respect to those of other galax- 
ies (Figure II. 4). Indeed, the usual distribution of magnitudes of old globular clusters 
follows a Gaussian law, which gives a log-normal shape for the mass function when as- 
suming a constant mass-to-light ratio (Harris 1991). In the Antennae however, the mass 
function seems to be best described by a power law M" of index (3 ~ —2 for the range 
10 4 M < M < 10 6 M (solid curves in Figure II.4). 

Theoretical studies suggested that mass-loss through two-body relaxation, tidal shocks 
and stellar evolution might be responsible for a dynamical transition of the mass function, 
from an universal power-law to a log-normal distribution, as observed in e.g. the Milky 
Way (Fall & Rees 1977; Vesperini 1998; Fall & Zhang 2001). This can be explained by the 
dissolution of the low mass associations (M < 10 4 M ) within the first 10 7-9 yr of their 
life, and a mild mass-loss of the more massive clusters, shifting to intermediate masses 
(10 4-5 M Q ). Ultimately (over ~ 10 10 yr), this leads to the mass function observed in the 
Milky Way for older objects. The effects responsible for the dissolution and mass-loss of 
the clusters can be either internal (violent relaxation) or external (via the tidal field), as 
presented in Section 3.3 (page 14). 

Fall, Chandar, & Whitmore (2005) studied the age of the clusters by comparing their 
magnitudes in different bands (U, B, V, I) with those from stellar population models 
(Bruzual & Chariot 2003). After applying a correction for completeness d and an age 

d Observational data are often biased by incompleteness, due to technical (e.g. limit in magnitude 
of the detector) or physical (like dust obscuration) reasons. Sometimes however, it is possible to correct 
for completeness (see e.g. Whitmore et al. 1999; Anders et al. 2007 for applications to the Antennae 
galaxies). 
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binning, they derived the feature-less age distribution dN/dr oc r 1 shown in green in 
Figure II. 5, with a median age of ~ 10 Myr. 

This can be interpreted by considering the birth and death rates of the clusters. The 
feature-less distribution is incompatible with discrete episodes of star formation, well 
separated in time, without any destruction. However, it is possible to assume a constant 
star formation rate 6 combined with a certain disruption rate. In order to get the r _1 
distribution, ~ 90% of the clusters must dissolve at each time decade. In other words, the 
younger bin is continuously filled by the constant formation rate. As time flows, 10% of its 
content is shifted to older age bins while the other 90% are destroyed. As a consequence, 
less and less clusters reach the oldest ages, which gives out the measured distribution. 
This rapid dissolution, often called "infant mortality" , has also been observed for clusters 

c This assumption is clearly an oversimplification, especially in the case of mergers, as noted by 
Whitmore et al. (20f0). The evolution of the (non-constant) SFR is discussed in Chapter III. 




Figure II. 4: Mass function of the Antennae star clusters, for two age bins. The 50% 
completeness limit is indicated by vertical dashed lines. The two solid lines represent 
power-laws ^(M) oc M 13 . The black dashed curve shows the log-normal mass function, 
for the clusters of the Milky Way, derived from a Gaussian distribution of magnitudes 
in the V band (< My >= -7.3 ; a(M v ) = 1.2, see Table 2 of Harris 1991) with a 
mass-to-light ratio M/Ly = 2. (Data points from Zhang & Fall 1999.) 
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in the Milky Way (Lada & Lada 2003) and the Small Magellanic Cloud (Chandar, Fall, & 
Whitmore 2006), over the first ~ 100 Myr. As mentioned before, this dissolution results 
from both internal and external effects, among which the tidal field of the host galaxy. 
That is why the age distribution of the star clusters is an important question that we 
will address in the following, from the standpoint of galaxy mergers. 



1.4 ... and tidal dwarf galaxies 

While star clusters are mainly located in the central region of the Antennae, it is worth 
noticing that other structures exist in the tidal tails. Schweizer (1978) and Mirabel et 
al. (1992) were the first to detect a pair of massive Hi clouds undergoing star formation, 
near the tip of the southern tail (see Figure II. 6). This has been interpreted as TDGs 
exhibiting properties similar to those of dwarf irregulars or blue compact dwarf galaxies 
(see e.g. Due & Mirabel 1999). 

These observations left open the question pertaining to the origin of these blue regions. 



T3 




10 6 10 7 10 8 10 9 

T [yr] 



Figure II. 5: Age distribution of the clusters in the Antennae, based on different mass 
criteria. The age binning is showed by horizontal bars. The green line represents a 
distribution dN/dr oc r _1 . (Data points from Fall, Chandar, & Whitmore 2005.) 
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Figure II. 6: Composite image of far-UV (blue) and near-UV (red) observations of the 
Antennae by GALEX (Martin et al. 2005), clearly showing the TDGs candidates near 
the tip of the southern tail. (Adaptation of Figure l.a of Hibbard et al. 2005.) 
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To address it, Hibbard et al. (2005) used UV data and optical images to study the ages of 
the stellar populations in these dwarf galaxies (Figure II. 7). They noted that a continuous 
formation of stars is very unlikely in these areas. The different densities of gas in the two 
tails have led to different star formation epochs which explains the spanning of relative 
age between the two tidal structures. 

SSP | In the southern tail, simple stellar population models reveal that most stars have 
formed before the tail itself (Hibbard et al. 2005). However, structures visible in UV yield 
younger ages than the dynamical age of the arm (estimated via numerical simulations), 
suggesting a recent episode of star formation. Therefore, the stars observed in the tails 




Figure II. 7: Intensity-weighted ages of the tidal regions derived from SSP thanks to 
(far-UV - near-UV) and (far-UV - V) data. A gradient is visible along the southern 
tail, from old populations close to the nucleus to younger near the tip. The pattern is 
more complex within the northern tail. The contour line represents a surface density of 
27.5 mag arcsec -2 in the near-UV band. The arrow on the color bar shows the dynamical 
age of the southern tail from the simulation of Barnes (1988, see the next Section). 
(Adaptation of Figure 2.c of Hibbard et al. 2005.) 
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are likely to be of two origins: (1) the material of the progenitor disks which has been 
expelled during the tidal interaction and (2) an in situ formation after the creation of the 
tails. As the bluest region, the TDGs are the youngest objects in the arm, and consist 
of two populations of young (4 — 30 Myr) and intermediate (80 — 100 Myr) ages. 

Both central and tail regions of the Antennae are pristine laboratories for the study 
of the formation of star clusters and TDGs. Clearly, the tidal field plays a major role in 
the evolution of these populations. Therefore, a 4D (space and time) description of the 
galaxies will bring important clues to interpret the observations. This is only feasible via 
numerical simulations. 



1.5 Previous simulations of the Antennae 



Since the first observations, the morphology of the Antennae raised many questions 
about their origin. In their precursor work, Toomre & Toomre (1972) made the first 
attempt to reproduce the galaxies with simulations (Figure II. 8, top panel). 

As mentioned in the Section 1.3 (page 4), these simulations are not self-consistent 
as they only compute the influence of two point-masses on pseudo-particles. However, 
this approach manages to reproduce the long tails and proves that they originate from a 
tidal interaction. The main problem of this model resides in its symmetry, visible in the 
extension of the tails. The authors also noted a too large separation between the nuclei 
and suggested to take orbital decay into account to reach a better configuration. Indeed, 
their point masses stay in their initial Keplerian orbit, as no phenomenon makes them 
lose energy. Such restricted simulations cannot model the transfer of energy through 
dynamical friction (see Section A. 2, page 198), as the pseudo-particles do not create any 
gravitational drag. 

The first self-consistent simulation of the Antennae came 16 years later, from Barnes 
(1988, see Figure II. 8, bottom panel). This model has been used to show the role of the 
dark matter halos, which act as "a sink for binding energy and angular momentum as 
the visible galaxies spiral together" , so that the nuclei are closer together. 

Mihos, Bothun, & Richstone (1993) performed the first simulation of the Antennae 
including hydrodynamics. With a gaseous component in addition to the stars and DM, 
the simulation shows a dissipative part, used to form new stars (see Figure II. 9). 
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Box II.2 

In hydrodynamic simulations, a "recipe" is needed to convert an element 
of gas into an element of stellar material. This prescription is often based on 
the empirical relation of Schmidt (1959), who established a link between the 
SFR per surface unit and the observed gas surface density S g : 

E* oc E" (11.19) 





Figure II. 8: Top: First simulation of the Antennae galaxies by Toomre & Toomre (1972). 
Filled or open symbols distinguish the origin of the particles. (Figure 23 of Toomre & 
Toomre 1972.) Bottom: First self-consistent simulation of the Antennae (Barnes 1988). 
Contrary to the restricted simulation of Toomre & Toomre (1972), this model accounts 
for the orbital decay and better reproduces the inner region of the merger. (Figure 2 of 
Barnes 1988.) 
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Figure II. 9: Top: First hydrodynamic simulation of the Antennae by Mihos, Bothun, 
& Richstone (1993). Middle: map of the star formation rate from the same simulation. 
(Adaptation of Figure 8 of Mihos, Bothun, & Richstone 1993.) Bottom: Recent SPH 
simulation by Karl et al. (2008). Old stars are green, new stars are red and the gas is 
shown in blue. (Adaptation of Figure 2 of Karl et al. 2008.) 
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This relation has been tuned by Kennicutt (1998), who determined 
n = 1.4 ±0.15. In the simulations, the Schmidt-Kennicutt law is converted to 
describe the SFR as a function of the gas volume density: 

p* (x Mp n g , (11.20) 

where M is the mass of the initial cloud (Mihos, Bothun, & Richstone 1993). 
The normalization factors are usually chosen to set the SFR of an isolated 
galaxy to ~ 1 M yr -1 . Other approaches (Katz 1992; Springel 2000) use the 
timescale for star formation, which is linked to the time needed by the gas 
cloud to collapse: 

P* « 777^' ( IL21 ) 
which gives a Schmidt law with n = 3/2. 



Such a method provides interesting clues about the location and age of newly formed 
structures. However, the physical details of the collapse of the molecular cloud and the 
other competing processes are not described, as noted by Barnes (2004). 

Many other models have been created, most of them based on the orbital parameters 
CPU | of Toomre & Toomre (1972). Today, some of the central process unit power is dedicated 
to the research of a better set of parameters, before performing the simulation itself (see 
e.g. Karl et al. 2008). In the next Section, we present our model of the Antennae and 
explain why such parameter studies are more than welcome! 



2 Numerical model 



In order to explore the tidal field of the Antennae, one needs a 3D representation of 
its potential. As discussed before, the gaseous component of galaxies does not strongly 
influence the gravitational potential, because of its relatively small mass. Therefore, stars 
and DM seem to be the only components required in the modeling of the gravitational 
field of the merger. As our simulations are collisionless, we use Nemo and gyrf alcON for 
the setup and the integration of the orbits (see Section 1.1.4, page 31 and Section 1.1.6, 
page 34). This Section details and justifies the choices made for the main parameters. 
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2.1 Individual galaxies 

During a merger, the tidal stripping mostly affects the dynamically cold components 
of the galaxies, i.e. their disk. Indeed, the observations of long tidal tails in the Antennae 
suggest that the progenitors are spiral galaxies, instead of ellipticals. Within the Nemo 
package, the tool magalie (Boily, Kroupa, & Penarrubia-Garrido 2001) creates a baryonic 
component made of an exponential disk and a Hernquist (1990) bulge embedded in an 
isothermal DM halo, in a similar way as Hernquist (1993) who presented a method to 
build multi-component models of galaxies. 

In the observations of the tidal tails, we noted an asymmetry between the two galaxies: 
the northern tail (associated with NGC 4039, the southern nucleus) shows a shorter 
extension than its counterpart. This suggests that the initial disk of NGC 4039 was 
more compact than the other one. However, the luminosity and size of the two nuclei are 
similar, which, in a first approximation, implies that they have similar masses. Hence, 
the radius of the disk of NGC 4039 has to be smaller than the one of NGC 4038. 

In our model, the halo component is made of DM only. However, it is clear that 
typical spiral galaxies also have a baryonic extension of a similar shape, called spheroid. 
Therefore, our bulge component from magalie is much more extended than a classical 
bulge, and can be mentally divided into an actual bulge and a spheroid of lower density. 
Thanks to this trick, all the baryonic matter is included in the disk and bulge components. 

The integration of the model of a single progenitor (i.e. in isolation) with gyrf alcON 
reveals that it is not initially virialized. Indeed, small scale relaxation effects occur for 
about a crossing time, while the overall morphology is not affected by the redistribution 
of the kinetic energy. To avoid any artifact, the isolated model is first left evolving to a 
stable virialized system, and then extracted to be used as initial configuration. 



2.2 Orbital parameters 

Once the parameters of the progenitors have been set, one can focus on the geometry 
of the merger. The aim here is to determine how the two galaxies interact with each 
other, in term of orbits, relative mass and size, inclination and so on. 

As presented above, the two progenitors have similar masses and thus, the Antennae 
are considered as a major merger (i.e. a pair of galaxies with a mass ratio greater than 
3:1). For simplicity, we assume that their rotation curves are also similar, which implies 
that no scaling (in mass, velocity or size) has to be done to the progenitors. We also 
keep the same parameters for the disks, bulges and halos (except the radius of the disk, 
as mentioned before). 

To determine the orbit of the galaxies, one has to set the initial conditions as differen- 
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tial position and velocity of the progenitors. Obviously, there are no strong constraints 
on the space of the orbital parameters because the only criteria are the observed veloci- 
ties and the morphology at present time, i.e. two quantities that depend on the line of 
sight and the initial inclination of the disks with respect to the orbital plane. That is, 
the only way to set the initial value of an orbital parameter is to assume that all the 
others are known, to integrate the model until it best matches the observations, and to 
iterate with another parameter. However, if one plays with, say, the orbital eccentricity, 
it is very likely that the inclination that would give the best match will change too. In 
other words, all the variables are inter-dependent and cannot be tuned sequentially, one 
after the other. 

Nevertheless, one has a few clues from observations and a basic knowledge of tidal 
interactions. For example, the large extension of both arms favors a prograde-prograde 
encounter (see Section 2.2, page 8). This constrains the inclination of the disks relatively 
to the orbital plane. Furthermore, the two nuclei are very close to each other, suggesting 
that they recently underwent or will soon undergo a peri center passage. Because the 
tidal tails had time to expand to their current positions, this passage is not the first 
one. After the first passage, the nuclei came back close to each other (as observed), 
meaning that the orbit is bound. However, one should not conclude that the initial 
eccentricity has to be smaller than unity. Indeed, in order to create the big tails, each 
progenitor underwent a close passage through the material of the other. Then, a fraction 
of the initial orbital energy has been transfered to the individual particles via dynamical 
friction (see an example in Figure A. 2, page 200). Therefore, even a Keplerian orbit with 
an eccentricity slightly greater than unity could lead to the bound orbit required by the 
Antennae data. 



2.3 Simulation of the merger 

It would be possible to explore the parameter space through a grid of models, but 
the number of combinations and the computational power needed to perform a single 
simulation forbid such a method f . One could also imagine to consider the observations 
as the initial state and go backward in time, up to a moment where the galaxies are 
isolated. Although possible in principle, such an approach would require a very accurate 
knowledge of quantities like the velocities perpendicular to the line of sight, or the mass 
distribution, which are not measurable via observations. 

Hence, the approach adopted here is an iterative method where each parameter is 
tuned manually until the simulations give an acceptable match to observations (mor- 
phology and kinematics). The quality of the match is determined by eye. Note that the 

f Several authors proposed to explore the parameter space with restricted simulations, that are 
obviously much faster to run than self-consistent ones (see e.g. Theis & Kohle 2001; Barnes & Hibbard 
2009). However, these simulations are still facing problems in the recovering of the orbital decay, an 
important property of mergers. 
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Table II. 1: Parameters of the Antennae model 



Parameter 


NGC 4038 


NGC 4039 


Numbers of particles 




^V"total 


7 x I0 5 


7 x 10 5 


Naisk 


2 x I0 5 


2 x 10 5 


-Abulgc 


I x I0 5 


I x 10 5 




4 x 10 5 


4 x 10 5 


Scalelengths 






^disk,vcrtical 


0.2 


0.2 


f disk 


1.0 


1.0 


'"bulge 


1.0 


1.0 


r halo 


7.0 


7.0 


Cut-off radii 






Cdisk,vertical 


2.0 


2.0 


Cdisk 


5.0 


3.0 


Cbulgc 


1.0 


1.0 


Chalo 


7.0 


7.0 


Masses 






"^disk 


1.0 


1.0 


"^bulge 


1.0 


1.0 


m halo 


5.0 


5.0 


Toomre parameters 






Q 


1.5 


1.5 


Disk inclinations 








60° 


-60° 


9y 


0° 


30° 


Initial coordinates 








(6.0,6.0) ( 


-6.0,-6.0) 


(vx,v y ) ( 


-0.5,-0.25) 


(0.5,0.25) 


Note: The parameters are given in numerical units (G = 1) with 


the conversion factors 4.4 


: kpc, 1.9 xlO 2 km s" 1 and3.6xl0 10 M 


for the length, velocity 


and mass respectively. 


The Keplerian 


equivalent orbit has an 


eccentricity of 0.96, a distance to peri- 


center of 1.31 (= 5.78 


kpc) and a semi-major 


axis of 17.74 



(= 78.06 kpc). 
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goal of this study is not to produce a better model for the Antennae than those from 
the literature, but simply to find a set of parameters which will, in the end give a good 
description of the tidal field. The final set of parameters, presented in Table II. 1, has 
been obtained after running ~ 100 low resolution attempts. Figure 11.10 and Figure 11.11 
display the morphology of the system, as the logarithm of the surface density, a good 
way to mimic the CCD images from observations (see Dubinski 2008 for an introduction 
to visualization of iV-body simulations). In fine, our model is similar to that of Barnes 
(1988), except that he used lighter King models for the bulges. The choice made here 
allows a finer description of the central region of the merger, closer to observations. 

The entire parameter space exploration and the simulations themselves are done us- 
ing numerical units. By comparing the final result with observations, it is possible to 
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Figure 11.10: Morphology of the Antennae model, shown in the orbital plane, every 
50 Myr. The shades of blue represent the logarithmic surface density, t = Myr corre- 
sponds to the first pericenter passage. 
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Figure 11.11: Same as Figure 11.10 but in the plane of the sky. The orbital plane is shown 
in red in the panel a. The arrowhead marks the position of the point (0, 0, +50 kpc) in 
the orbital plane (i.e. foreground), which is now behind the plane of the sky (the arrow 
is pointing away from the reader). The best match with the observations is obtained for 
t = 300 Myr (panel h, see also Figure 11.12). 
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determine the conversion factors. Observational data constrain well the size of an object 
and its velocity along the line of sight. The mass and time however, are not directly 
measured. That is why the non-dimensionalization used here is based on the spatial 
extension of the tidal tails and the peak in the velocity field measured by Hibbard et al. 
(2001). With these factors and the relations given in Section 1.1.2 (page 26), one can 
retrieve the equivalent factors for the mass and time, and deduce the physical resolution 
of the simulation. In this case, each particle represents ~ 2 x 10 5 M Q and a snapshot 
corresponds to 2.5 Myr. (Hereafter, the origin of time t = corresponds to the first 
pericenter passage.) The smoothing length used covers a radius of 40 pc. These num- 
bers match well the physical size and mass of a typical star cluster, which means that a 
particle in this simulation represents a single cluster (whose internal dynamics are out of 
reach) . 

The best match with observations is reached 300 Myr after the first pericenter 
(Figure Il.ll.h and Figure 11.12). Despite the too long northern tail, the model dis- 
plays a central region showing most of the structures of visible in e.g. HST images. The 
velocity field (along the line of sight) is shown in Figure 11.13. The general structures are 
well reproduced, notably the gradient along the southern tail, up to the tip which shows 
a "folding" pattern. We also remark that each nucleus yields a v z component of opposite 
sign to those of the related tail, as previously noted by the observations (see Figure II. 3). 
On a smaller scale, the simulation fails to show the internal rotation of the nuclei. A 
multicomponent distribution is visible in the central region, showing the complexity of 
the system that is currently undergoing an efficient mixing of the matter. 

The good match with observations for both the morphology and the kinematics gives 
us confidence in this model for the next steps of our study. Before that, we need to 
characterize the simulation with more details, in particular its evolution in time. 

Figure 11.14 shows the orbits of the progenitors in the orbital plane, corresponding to 
the trajectories of the center of mass of the most bound particles in each galaxy. Indeed, 
by considering the center of mass of all the particles, one would take into account the tidal 
tails, which would lead to an offset in position while the nuclei merge (see an example 
on Figure A. 3, page 201). This effect disappears when selecting the particles that stay 
in the nuclei throughout the simulation^ As a complement, Figure 11.15 displays the 
evolution of the three-dimensional distance between the two centers of mass. The two 
pericenter passages are clearly visible as local minima for the distance. This plot tells 
us that the current configuration (t = 300 Myr) is reached just after the second passage. 
One can also remark that the dynamical friction forbids a third separation of the two 
nuclei, which rapidly merge. The evolution of the distance between galaxies in mergers 
is of paramount importance when considering the SFR or the effect of the tides, as we 
will see below. 



g The center of density or potential give good results too, even when considering all the particles. 
However these quantities are determined at high computational cost. 
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Figure 11.12: Morphology of the Antennae model for t = 300 Myr, at large scale (left), 
and in the central region (right). The southern tail matches well the observations, al- 
though the northern one is too extended. In the central region, the origin of the tails is 
clearly visible, as well as the main features: the two nuclei, the overlap region and the 
arc shape in the northern part. 



3 Tidal field of the progenitors 



The main topic of this thesis is to quantify the effect of the tides on the baryonic 
matter. For this reason, we concentrate our exploration on the stellar particles, i.e. the 
disks. The results presented in the following make intensive use of the method to compute 
the tidal tensor at the position of the particles of the disks (recall Section 1.2, page 38). 
With the resolution adopted here, the cubes used for the finite difference scheme have a 
size 5 ~ 220 pc, i.e. well beyond the half-mass radius of typical clusters. 

Before studying the tidal field of the merger, it is important to get familiar with those 
of the progenitors. Indeed, if every object creates its own tidal field, it is the gathering of 
two of them that will build an interesting pattern, responsible for e.g. the formation of 
the tidal tails. That is why, in our approach, it is important to distinguish the intrinsic 
tidal field from the one due to external objects (a second galaxy in this case). It is also 
a good opportunity to become familiar with the tools and diagnostics we will use to 
study the statistics of the tidal field of much more complex systems. In this Section, we 
consider the model of NGC 4038 in isolation. Note that, despite a small difference in the 
radial parameter of the disk, NGC 4039 yields comparable results. 
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Figure 11.13: The velocity field along the line of sight of the Antennae model (left) and 
a zoom-in on the central region (right), at present time. The dashed square represents 
approximately the size of the right panel of Figure II. 3. The kinematics of the tails are 
well reproduced by the simulation but the internal rotation of the nuclei is not clearly 
visible. Note that the peak velocity matches de facto the observations because it has 
been chosen for the non-dimensionalization of the simulation. Therefore the quality of 
the model is estimated from the gradients along the tails, and in the central region. 

3.1 Compressive nucleus 

Figure 11.16 (top panel) shows the radial profile of the potential <fi of one of the progen- 
itors of the Antennae model, taken in isolation. The associated mathematical functional 
form shows a change of curvature from convex to concave at a radius r c ~ 1.3 kpc (verti- 
cal red line). As the second derivative of the potential, the tidal tensor T is the negative 
of its Hessian matrix, so a change of curvature in implies a change of sign in T. There- 
fore, in its diagonal form, the dominant term of T goes from negative to positive with 
increasing radius. The other eigenvalues being negative everywhere, the tides are fully 
compressive (see Section 2.4, page 11) in the inner w 1.3 kpc region. On the bottom 
panel of Figure 11.16, one reads that the corresponding mass represents ~ 2.5% of the 
total mass of the disk (i.e. ~ 9.0 x 10 8 M ). 



3.2 Three classes of orbits 

These figures depend on the shape of potential of the progenitor, and thus should not 
vary with time. However, it is important to understand that only the statistics of the 
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Figure 11.14: Trajectories of the centers of mass of the most bound particles of the 
two progenitors, again in the orbital plane. The pericenter passages are marked with • 
(t — Myr, Figure Il.lO.b) and + (t — 285 Myr). The configuration closest to observa- 
tions is shown with * (t = 300 Myr, Figure Il.lO.h). 
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compressive modes are constant but the evolution of the tidal field of individual particles 
is not. In a disk, the mean motion of the particles is the rotation. But radial and vertical 
components also exist, voiding the orbit of a given particle to be a perfect circle. In this 
Section, we detail the impact of the orbital parameters of the particles into the statistics 
of the tidal field, in term of duration of the compressive mode. 

In projection on the plane of the disk, the possible cases for orbits are illustrated in 
Figure 11.17. The particle on orbit A has an energy and angular momentum high enough 
to stay permanently outside of the compressive region. On the orbit B, the particle shows 
a similar level of energy, but a lower angular momentum. That is why it has a small 
pericenter distance and thus, enters and leaves the compressive area. Finally, the orbit 
C has a low energy level and remains in the compressive region. 

Therefore at a given time t, the fraction of the particles in compressive mode (i.e. closer 
than ~ 1.3 kpc to the center of the galaxy) are made of large binding energy (class C) 
and lower binding energy but low angular momentum (class B) particles. At t + dt, all the 
particles of the class C are still in this compressive region but some of the class B have left 




Figure 11.15: Evolution of the distance between the two progenitor galaxies. Again, the 
pericenter passages and the current configuration are marked with •, + and respec- 
tively. The merger phase is clearly visible after t ~ 350 Myr, when the distance shows a 
damped oscillatory evolution. 
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or entered it. Statistically, the total mass fraction in compressive mode remains ~ 2.5% 
but the particles which build it have partially changed. Figure 11.18 plots the total 
energy E and the angular momentum L of the particles of the disk, as in Hemsendorf, 
Sigurdsson, & Spurzem (2002). One reads the minimum energy E c w —6.9 x 10 39 erg 
and the maximum angular momentum L c w 5.0 x 10 53 erg s at the boundary of the 
compressive region (r c ~ 1.3 kpc). With these values, it is possible to plot the diagnostic 
Figure 11.19, and to distinguish the three classes of orbits. 

In fine, when considering the timescale of the compressive mode, one has to take into 
account the contribution of the three classes: class A yields a null period in compressive 
mode, class B shows a characteristic time r, depending on the details of the orbits of 
its particles, and class C has an arbitrarily long time in compressive mode. To estimate 
r, we compute the free-fall time of a particle through the compressive region and get 
r ?y 37 Myr. We note that this value is well approximated by the period of a quarter of 




5 10 15 20 

r [kpc] 



Figure 11.16: Top: Gravitational potential of the disk particles of one progenitor galaxy. 
A change in the curvature at r c f=s 1.3 kpc is marked by a red line: using mathematics 
terminology, the potential <f>(r) is a convex function for r < r c , while it is concave 
elsewhere. In other words, the second derivative is positive in the inner region, which 
corresponds to a compressive tidal field. Bottom: Normalized integrated mass of the 
disk. The radius r c w 1.3 kpc encloses ~ 2.5% of the total mass. 



77 



Chap. II 



The Antennae galaxies 




Figure 11.17: The three possible cases of orbits of particles in a disk, with respect to 
the compressive radius (dashed black circle). (The eccentricities have been exaggerated 
for ease of illustration.) Orbit A (blue) is entirely outside the radius of the compressive 
region while orbit B (green) enters and leaves it. Orbit C is enclosed by the compressive 
radius. The x sign marks the center of the disk. Note that the vertical components of 
the orbits are not taken into account explicitly, for simplicity. 



circular orbit at the radius r c : 



It is important to make the difference between the time needed by a particle 
left at r = r c without an initial velocity to reach r = 0, and the time the sphere 
of radius r c will take to collapse if its kinetic energy is instantaneously set to 
zero. In the first case the mass of the inner region (source of gravity) depends 
on the position of the particle, and thus evolves with time, while in the later 
case, this mass is kept constant. For a homogeneous sphere of density p , one 
can show that the first timescale is 



T ■ 

± r.i 



4 4^GM(r < r c ) 
« 37 Myr. 



(11.22) 



■ Box II.3 




(11.23) 



and the second (more often used in the literature) is 
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(11.24) 



In our case, we are interested in characterizing the motion of a particle inside 
the compressive region, thus the first timescale. 



The real time spent in compressive mode must be shorter than the time needed by a 
particle with no initial velocity to get from r c to the center and then to r c again (enter, 
cross and exit the compressive zone), i.e. approximately twice our estimate of the free-fall 
time: ~ 75 Myr. 

This can be interpreted by considering the circular orbit of a particle situated at r c , 
plus an epicyclic variation which makes it oscillate around its mean radius. That is, we 




Figure 11.18: Total energy (top) and angular momentum (bottom) of the disk particles. 
The vertical red line marks the boundary r c of the compressive region. Horizontal lines 
show the level of energy E c and angular momentum L c at this boundary, used in the 
following diagnostic. The green curve represents the theoretical angular momentum of 
circular orbits. 
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Figure 11.19: The three classes of orbits A (blue), B (green) and C (red) are distinguish- 
able in the energy- angular momentum plane. (Note that only the particles within the 
central 5 kpc are displayed.) The value of energy E c and angular momentum L c at the 
boundary of the compressive region are marked with dashed lines. The class A particles 
have E > E c and L > L c and thus never enter into the compressive region. The B 
class gathers particles with a high energy level (E > E c ) but a low angular momentum 
(L < L c ) which allows them to enter the central region for a finite amount of time. The 
C particles have a low energy (E < E c ) and low angular momentum (L < L c ) and hence, 
remain in the compressive volume. For each class, the percentage indicates the fraction 
of the total disk mass. 
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characterize the fraction of the orbital motion during which the particle goes inside the 
compressive region, over the total orbital period. 

In this simple approach, we have neglected the two-body relaxation, the orbital decay 
and possible resonance phenomena. Including these second order effects would modify 
our conclusions for some particles but not for the overall statistics. This assumption 
has been confirmed by the relative stability of the previous results, when applying the 
same methods to different snapshots (and thus different epochs). Figure 11.20 shows the 
evolution of the mass fraction in compressive mode over 10 Gyr, i.e. much longer than 
the dynamical times of merging systems. A systematic increase (0.05% Gyr" 1 ) is visible 
however, and can be attributed to secular evolution of the galaxy (possibly a slow transfer 
of angular momentum through a weak bar mode). In the following, this weak effect has 
been neglected. 




Figure 11.20: Long term evolution of the mass fraction in compressive mode. (Here again 
the time resolution is 2.5 Myr.) Over 10 Gyr, the fraction changes by a few 0.1% only, 
from its initial value. 
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3.3 Statistical timescale of compressive modes 

One seeks a statistical description of the period spent in compressive mode to quan- 
tify the idea introduced above. In Renaud et al. (2008), we defined two quantities to 
distinguish between the two characteristics of the time evolution of the tidal field: the 
longest uninterrupted sequence (lus) and the total time (tt). For a given particle, lus is 
the duration of the longest continuous episode in compressive mode, while tt is the sum 
of all compressive episodes, i.e. the total time during which the particle stands in a 
compressive tidal field. For example, consider a particle in extensive mode for 80 Myr, 
going to compressive mode for 40 Myr, extensive again for 30 Myr and finally compres- 
sive during 50 Myr. Its lus is 50 Myr while its tt equals 40 + 50 = 90 Myr. In any case, 
one has lus < tt, the equality being reached when the particle experiences only one or 
zero compressive event. Therefore, the three classes of orbits yield different statistical 
timescales: class A shows lus = tt = 0, class B has < lus < tt and class C yields the 
largest values of lus (which also equals tt). 

Figure 11.21 plots the distributions of lus and tt, in term of mass fraction*: for a 
given time interval At, the curve gives the mass fraction of the disk that yields a lus 
(or tt) longer than or equal to At. The plot is normalized by considering only the 
non- constantly compressive particles (i.e. classes A and B). Therefore, the curves on 
these plots are always decreasing with a maximum value of 100% for At = (since any 
compressive mode lasts more than Myr). 

In Figure 11.21, the lus and tt distributions are equal for the period of At < 2.5 Myr 
(i.e. one snapshot, the time resolution) and ~ 92% of the mass. This matches well our 
estimates for the class A. The remaining ~ 8% are the class B particles, with a lus longer 
than one snapshot. The long periods reflect the fluctuations previously mentioned (e.g. 
two-body relaxation) and thus corresponds only to a small number of particles. 

Both distributions can be fitted over the interval ]0, 200] Myr by double exponential 
laws of the form 



where fi and f2 are constants, and T\ and T2 are characteristic timescales. Table II. 2 lists 
the best fit parameters and the residual error a, as shown with dashed line in Figure 11.21. 
These double exponential laws suggest two regimes defined by the time intervals ATi and 



h The reader might take the time to get familiar with this kind of figure, as it will be often used in 
the next Sections and Chapters. 




(11.25) 
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AT 2 during which one exponential dominates the other J . For the lus distribution, the 
transition occurs at At ~ 27 Myr, where the distribution shows a knee in the log-normal 
plane. 

Note that these results depend on the duration of integration. Indeed, when integrat- 
ing over a long time, one allows the particles to stand more often in a compressive mode. 
For quiescent systems, like our isolated progenitor, the lus is set by quantities that are 
roughly conserved over time (E and L). With no perturbation, the particles remain on 
similar orbits and thus have the same lus for each revolution. The tt however necessarily 
increases with time and leads to a flatter distribution. Figure 11.22 illustrates this point 
by showing the lus and tt distributions for different integration periods. The lus varies for 



1 In Renaud et al. (2008) and Renaud et al. (2009), the distributions have been split into short 
and long term regimes, and each has been fitted separately. In the present document, the distributions 
are fitted with Equation 11.25 over the entire range ]0, 200] Myr. This avoids the approximation in the 
definition of the transition zone between ATi and AT2, and generally gives better fits. Therefore, the 
parameters presented here differ from those in our previous publications. The overall distributions are 
similarly reproduced however, and the conclusions in all documents are not affected by this change. 
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Figure 11.21: lus and tt distributions of compressive mode, integrated for 1.1 Gyr. The ar- 
row indicates the time interval where the class B particles become important. The dashed 
lines corresponds to the double exponential fit to the numerical results (see Equation 11.25 
and Table II.2). 
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Table II. 2: Parameters of the lus and tt distributions 



Distribution 


fi 


n [Myr] 


h 


T 2 [Myr] 


a 


lus 


10.1 


6.2 


0.2 


69.9 


0.16 


tt 


5.0 


7.6 


3.5 


216.1 


0.15 



no more than 1% while differences for tt reach one order of magnitude for At = 200 Myr. 
In any case however, the distributions are still well fitted with double exponential laws 
and the change of regime between AT\ and AT2 remains at approximately the same At. 




Figure 11.22: lus and tt distributions of compressive mode, for several integration times 
(0.5 Gyr, 1.0 Gyr and 2.0 Gyr). As expected, the long integration periods yield a higher 
integrated fraction, as one gives more opportunities to the particles to go into the com- 
pressive region. 
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* The Id letters refer to the panel identifiers in 
Figure 11.10, Figure 11.11 and Figure 11.25, as well 
as time markers in Figure 11.23. 



4 Tidal field of the merger 

In the previous Section, we have characterized the properties of the tidal field of the 
progenitor galaxies of our Antennae model. It is now possible to apply the same method 
to the merger itself and quantify its effect on the tidal field, with respect to isolated 
galaxies. In a first approach, we focus on a binary test for the nature of the tides: 
compressive or extensive. To do so, the sign of the maximum eigenvalue of the tidal 
tensor for all the particles of both disks is examined. 



4.1 Mass fraction in compressive mode 

For an isolated progenitor, the mass fraction in compressive mode depends on the 
number of particles lying in the central region, i.e. ~ 2.5%. We have seen that this 
number is statistically relevant and remains stable for an arbitrarily long time. During 
the course of the merger however, the fast evolving shape of the global potential is 
expected to influence this value by creating local, off-nuclear cores and thus additional 
compressive regions. 

Figure 11.23 shows the evolution of the mass fraction in compressive mode, from the 
isolation stage to the merger phase. The relation between major interaction events and 
the increase of this fraction is clearly shown in this plot as peaks (see also Table II. 3). 
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Both passages lead to a rapicV increase of the mass fraction in compressive mode, by 
a factor rs 5 compared to the level in isolation. The decreasing rate then depends on 
the configuration of the merger. When the galaxies fly away from each other (steps d, 
e, f and g), the fraction is slowly and smoothly decreasing, while it changes much more 
quickly during the merger phase (step i) . This last phase is well visible as a fairly regular 
damped oscillation. We note that the fraction is comparable for both isolation regimes: 
before interaction (t < —20 Myr) and after the merger (t > 450 Myr). 

As for an isolated progenitor, the final fraction (~ 3%) remains stable for a very long 
time. It maintains however a mild oscillation (mostly covered by the noise) with the same 
period than the oscillatory merger phase (~ 20 Myr), well visible in the frequency do- 
PSD | main: Figure 11.24 shows the power spectral density of the mass fraction in compressive 
mode. 



J The rise begins a few Myr before the pericenter passage itself. A compressive event appears just 
before the distance reaches its minimum, because of the extended size of the disks. The shift in time 
between the beginning of the compressive peak and the maximum corresponds to the time needed for a 
galaxy to fly from the edge to the center of the other. 




Figure 11.23: Evolution of the mass fraction in compressive mode through the course of 
the merger of the Antennae galaxies. The two blue dashed lines indicate the first and 
second pericenter passages. 
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4.2 Spatial distribution of the compressive tidal mode 

In isolation, the tidal field shows a compressive mode in the central region of the 
galaxy. In the course of a merger, the tidal forces responsible for the creation of the large 
scale structures like bridges and tails are likely to disturb the character of the tidal field 
on similar scales and build new compressive regions. 

Figure 11.25 displays the morphology of the merger and marks particles standing in 
compressive mode with red dots. The panel (a) clearly shows the intrinsic compressive 
region discussed above (r c « 1.3 kpc), as the progenitors can still be considered isolated. 
During the first pericenter passage (b), the compressive regions expand to the local cores 
in the potential, where the two disks overlap. The tails that grow while the progenitors 
fly away from each other (panels c, d, e and f) contain several compressive regions, well 
separated from the nuclei of their host galaxy. The largest of those is situated near the 
tip of both tails, i.e. close to the position of the observed TDGs (recall Figure II. 6). 




Figure 11.24: PSD of the mass fraction in compressive mode (red) after the merger phase 
(t G [475,4725] Myr). The main peak at 20 Myr (surrounded by a forest of smaller 
yet significant ones) corresponds to the oscillation of the two nuclei, already visible in 
previous figures. Secondary peaks (52, 75, 125, 145 Myr) are also visible in the PSD of 
the 10% Lagrange radius (blue), which indicates a "breathing" of the merger remnant. 
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Figure 11.25: Spatial distribution of the particles in compressive tidal mode in the Anten- 
nae model, in the plane of the sky. The particles in compressive mode are marked with 
red dots. On each row, the black line represents 10 kpc. Again for the sake of clarity, 
only one out of 40 red dots is displayed. (Figure 3 of Renaud et al. 2009.) 
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Note that the tidal bridges also display compressive areas, but significantly less than 
the tails. A transient spiral pattern already existing in the morphology of the merger 
is partially reproduced in the map of the compressive mode for t « 200 Myr (see the 
inverted S-shape near the southern nucleus in panel f). 

During the second passage (panels g and h), the overlap of the two disks rises a 
large compressive region (see Figure 11.26). Because of the loss of their orbital angular 
momentum, the two galaxies spiral in on each other, leading to an oscillatory phase visible 
in Figure 11.15. The same happens for the compressive regions. At the end of the merger 
phase, we note the presence of two zones, visible on Figure II.25.i, expanding radially in 
a similar fashion than a wave. Finally, all the off-center compressive areas vanish within 
a few Myr, but the newly formed central nucleus remains tidally compressive. 




Figure 11.26: Map of the central compressive regions for t = 270 Myr, i.e. 10 Myr before 
the second pericenter passage. The spiral arms formed during the first passage are clearly 
visible around the two nuclei, as pseudo loops in this projection. 
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It appears that, during the course of the merger, compressive regions match the loca- 
tions of over-densities (overlap, nuclei) but also intermediate and low densities (bridges, 
tails). Figure 11.27 displays the ratio between the column density of compressive particles 
and the one of all particles, for NGC 4038 at t — 150 Myr. The gradients visible in this 
plot denotes that the compressive mode is independent of the density in the tidal struc- 
tures but shows a correlation with the densest regions in the nuclei (a similar pattern 
also exist for the northern galaxy). 

Figure 11.28 does not only shows the compressive regions but the entire tidal field, 
computed in a regular grid, for the central part of the Antennae. It also overlays the 
positions the young cluster candidates observed by Mengel et al. (2005). Most of these 
clusters lie in the red, compressive areas of the map, i.e. in the overlap region and the 
western-loop. Even when taking the projection effect into account, the strong correlation 




Figure 11.27: Ratio of the column density of the compressive particles S c to the total 
column density £ of the southern galaxy in the plane of the sky at t = 150 Myr (see the 
equivalent morphology in the thumbnail). The column densities are computed on a grid 
made of 300 pc-wide squared cells. In the nucleus, most of the particles are compressive 
but the fraction is much smaller in the tidal tail. Note that a ratio of unity is reach for 
the volume density in the nucleus, but not for the column density, because of foreground 
and background particles that are not in compressive mode. 
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visible between the young clusters and the compressive regions points out the role of the 
tidal field in the early phases of the life of a star cluster. Still the question of the timescale 
of these modes persists. 



4.3 Age distribution of the compressive tides 

One may argue that the good match between the positions of the observed young 
clusters and the existence of compressive mode might be a coincidence. A comparison 
between the duration of both phenomena (compressive tidal field and star formation) 
would be another point to confirm or not the supposed link between the two. To do so, 




Figure 11.28: Tidal field of the central part of the Antennae galaxies for t = 300 Myr. 
The blue to red colors label extensive to progressively more compressive mode. White 
corresponds to a neutral tidal field. The black dots indicate the positions of the young 
cluster candidates identified by Mengel et al. (2005, see their Figure 1). The map covers 
a ~ 12 kpc wide square, parallel to the plane of the sky, « 2.2 kpc backward from the 
plane containing the global center of mass. The nuclei of the two galaxies are marked 
with yellow crosses. Arrows roughly indicate the two major off-nuclear regions: the 
overlap and the western-loop. (Adaptation of Figure 2 of Renaud et al. 2008.) 
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Table II. 4: Parameters of the lus distributions 



Region 



fi 



n [Myr] h r 2 [Myr] 



Isolation 
Antennae 
Southern nucleus 
Southern arm 
Bridges 



10.1 
137.2 
32.1 
12.3 
21.3 



6.2 0.2 69.9 
6.8 0.2 119.3 

7.3 0.1 120.3 
7.4 



5.6 



we compare the time during which this mode is operative, to the typical timescales for 
the formation of the star clusters, i.e. a few Myr. We do this statistically, thanks to the 
lus and tt previously introduced. 

Figure 11.29 compares the lus distribution of the merger with the one of an iso- 
lated progenitor. The interaction clearly induces more short-term compressive events 
(At < 50 Myr) than for the progenitors taken in isolation. Here again, the lus dis- 
tribution can be fitted with a double exponential law (recall Equation 11.25) with the 
parameters given in Table II. 4. The timescale for the first period (rj) remains very simi- 
lar to the isolated case (only the mass fraction changes) , however r 2 is significantly higher 
for the merger. Note that the position of the knee marking the transition between AT\ 
and AT 2 is shifted from pa 27 Myr to pa 50 Myr. All together, these points underline that 
the interaction drives many short-term compressive events but also a few much longer 
episodes than in a quiescent disk. 

The underlying reason for these differences resides in the interaction between the two 
galaxies, and thus depends on the tidal structures. To investigate this idea, different 
regions of the merger have been explored separately. Clearly, the double exponential 
distribution is due to the nuclei of the galaxies which yield long compressive episodes. 
The short-term features on the other hand come mainly from the transient tidal structures 
like the tidal tails and bridges, which do not show any compressive episodes over AT 2 . 
The shorter timescales T\ however, are of the same order of magnitude everywhere. 

It is worth noticing that the periods reported here (t\ ~ 10 7 yr, knee at At pa 50 Myr) 
match well the timescales of star cluster formation. These numbers even suggest that the 
compressive mode could be active as long as the early phase of the life of a young cluster, 
before the tidal field switches back to the more usual extensive mode. The consequences 
of this change and their links with the very survival of these stellar associations will be 
discussed further. 
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Compressive mode interval At [Myr] 



Figure 11.29: lus {top) and tt {bottom) distributions in the Antennae (solid black line), 
compared to a progenitor in isolation (dashed black, as in Figure 11.21). Several regions 
of the Antennae have also been extracted, as shows in the thumbnail (for t = 150 Myr). 
For these subsets, the integrated fraction is normalized to the entire merger. Indeed, the 
value read for At = is the mass fraction corresponding to the selected region: 19% for 
the southern nucleus (NGC 4038), 13% for the southern arm and 14% for the bridges. 
Note that the part of the arm closest to the nucleus has not been selected because it 
rapidly falls back into the central region and thus, biases the statistics. In all cases, the 
integration has been performed for t G [—100, 1000] Myr. 
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4.4 Realizations of the statistics 

We have introduced the concepts of lus and tt distribution to get a statistical de- 
scription of the duration of the compressive tidal mode, applied to each mass element 
of our simulations. However, any statistics is a collection made of particular cases and 
it is important to keep in mind that the actual behavior of the tidal field for a given 
particle may not be well represented by the statistics presented above. In the following, 
we extract some typical cases for the evolution of the tidal field, in order to complete the 
picture we have built with statistic tools. 

Figure 11.30 displays a sample of seven orbits, with dots marking the compressive 
events along the orbits. Particles labeled A and B stay in the vicinity of the nucleus of 
NGC 4038. Particle A even remains in the intrinsic compressive region of the galaxy and 
thus, shows an almost uninterrupted compressive episode. By contrast, B only yields 
a relatively short compressive event, during the first pericenter passage. C goes to the 
bridges before falling back into the nucleus. Similarly, D falls back during the merger 
phase, an episode clearly denoted by a compressive sequence after the second passage. 
On a larger scale, E is ejected to the tidal arm by the first interaction but is captured 
back by the nucleus before the second. F and G remain in the tidal tails, with G being 
part of the TDGs region spotted in Section II. 4. 2. 

Such a plot shows the tidal history of individual particles or, in other words, some 
realizations of the statistics presented before. The major families of orbits are well 
summarized in this Figure which allows us to get a broader picture of how the tidal field 
evolves. In addition, this kind of approach is the starting point of the simulation of star 
clusters, using the tidal field derived here as boundary conditions, as presented in the 
next Chapter. 
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Figure 11.30: Sample of orbits (in the orbital plane) of seven particles in the central 
region (left) and at large scale (right). A compressive event for a given particle is shown 
by a dot on its orbit. The initial position of each particle is labeled by a letter. The 
pericenter passages of the progenitors are marked with crosses (+ and x, respectively) 
on every orbit. 
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Chapter hi 

Tides and star formation 



Overview ^ 

In this Chapter, the results on the tidal field of the An- 
tennae galaxies are applied to the context of the forma- 
tion of star clusters. We highlight the links between the 
scalelengths of star formation and those of the galactic 
tides, in which the clusters are embedded. Some of this 
work has been presented in Renaud et al. (2008) and 
Renaud et al. (2009). 

J 



1 Timing and energy of a compressive mode 

In the previous Chapter, we have seen that the compressive mode of the tidal field 
shows properties similar to those of the formation of star clusters and TDGs in the 
Antennae galaxies. First, the location of the compressive regions in our simulation cor- 
responds surprisingly well with the sites of numerous candidate young clusters. Second, 
the duration of these modes is of the same order as the timescale of star formation. 
These correlations are a strong hint that the tidal field, and especially its compressive 
modes, play a significant role in the formation of star clusters. It is therefore of prime 
importance to compare quantitatively the tides with other mechanisms involved in the 
formation of clusters. 



1.1 The compressive mode as a destructive effect 

As already suggested, the tidal field plays a role in the formation of the star clusters. 
It may also participate in their destruction either through: 
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• the expansion of the external layers triggered by gas mass-loss and enhanced by 
extensive tidal mode, 

• tidal shocks, as introduced by Spitzer (1958, see also Spitzer 1987, Kundic & Os- 
triker 1995 and Gnedin, Hernquist, & Ostriker 1999). 



Box IH.l 

In the context of star clusters, a tidal shock was first described by Spitzer 
(1958) as the dynamical heating of a cluster going through a galactic disk. 
Let's consider the density profile of a radially infinite disk, i.e. only depending 
on the height. A cluster centered on Z c undergoes the acceleration g(Z c ) from 
the disk. A star at the height z with respect to the center of the cluster [ergo 
at Z = z + Z c with respect to the disk) has the acceleration g(Z). In the 
reference frame of the cluster, its acceleration is therefore: 



dg_ 
dZ 



« z 



{111.1) 



The velocity of the cluster (~ 100 km s" 1 ) and the thinness of the disk 
(~ 100 pc) yield a crossing time of ~ 1 Myr, i.e. much shorter than the 
crossing time of the stars within the cluster. That is, it is possible to neglect 
the proper motion of the stars during the tidal shock, which is the framework 
of the impulsive approximation. Therefore, z remains constant and dZ can be 
replaced Vdt (where V is the constant orbital velocity) which gives 

Av z = ^Ag 

= -2<? m ^, (in.2) 

where g m is the value of the acceleration taken at large heights (i.e. where the 
density is low). The velocity after the shocks becomes 

v' = v + Av z , (III.3) 

and the new kinetic energy per unit mass is 

E' = X - (v 2 + 2v ■ Av z + Av z 2 ) 
AE = v • Av z + ^Av z 2 . (III.4) 
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The first term can be rewritten as v z Av z and yields a null average because 
of the symmetry in the distribution of v z . As a result, the average change in 
energy becomes 

(AE) = 2g 2 m ^ > 0, (III.5) 

which translates into a gain in energy for the stars, thus a dynamical heating 
of the cluster. This approach also applies when a cluster passes through a 
spiral arm, as noted by Gieles, Athanassoula, & Portegies Zwart (2007). 



Note that when a cluster is tidally shocked (as in the formalism of Spitzer 1987), the 
tidal forces only exist along one dimension and are actually compressive. However, the 
gain in energy of the cluster members does not depends on the sign of the tidal field and 
thus, the compressive modes have a destructive effect in this impulsive fashion. In the 
adiabatic regime however, the compressive modes act as a protective cocoon, by contrast 
with the extensive modes which enhance the dissolution (see Section 2.2 and Section 3 
below) . 

1.2 Age of the compressive events 

The statistical study of the lus- and tt distributions for the isolated progenitors and 
the Antennae is key to understand and describe the evolution of the tidal field. However, 
as noted in the previous Chapter, the actual history of an individual particle is a much 
more complex realization of these regimes: a particle can enter the compressive regime at 
the time of, say, a pericenter passage and remain in it for a certain period of time. After 
that, it may enter an extensive region for e.g. several xlO 7 yr, before immersing itself in 
a compressive mode again, and so on (recall the examples given in Figure 11.30, page 95). 
As a result, the interpretation of lus- and tt profiles independently is not sufficient to 
determine the net effect of the tidal field on a given particle along its orbit. As shown 
below, a combination of the two distributions provides a better statistics of the total age 
of the compressive mode that mimics the repeated switch from extensive to compressive 
modes. 

The age distributions of both lus and tt have been well-fitted with double exponential 
laws, with short (n) and long (r 2 ) timescales, corresponding to an early regime A7\ 
and a later phase AT 2 . A plausible interpretation is that the compressive history is 
first continuous during the short interval A2\. Then, when the tidal field evolves more 
significantly during the longer time lapse AT 2 , the tides go through a succession of discrete 
compressive episodes, described by the tt statistics. Hence, it would seem likely that the 
compressive event is continuous over the short timescale of lus (i.e. ti US) i), and that it 
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breaks up into several episodes over the longer timescale of tt (i.e. T tt 2 ). Consequently, 
the proxy for the actual time spent in compressive mode is a combination 3 of the lus 
distribution over the interval AT\, and of the tt distribution over the interval AT 2 : 



/ ie xp +/ 2 exp ). (III.6) 

TWl/ V r tt,2, 



Figure III. 1 displays the lus- and tt distributions plus the new combination in the 
log-log plane. The latter can be fitted with a power law dN/dt oc t~ L02 , over the range 
[5, 100] Myr (with a residual ~ 10~ 2 ). We note that the power law does not match our 
combination for ages t > 300 Myr, which corresponds to the time lapse between the two 

a This differs from the results presented in Renaud et al. (2009), where we have explicitly split the 
lus and tt distributions into a short (]0,50] Myr) and a long ([50, 100] Myr) timescale. This has given 
two single exponential fits per distribution. The approach followed here does not make this separation 
as our distributions are directly fitted by double exponential laws over the entire time interval. The 
separation between the interval is no longer arbitrary and thus, the quality of the fits is better with this 
new method. 




Figure III. 1 : Age distribution of the compressive modes, as a combination of the lus and 
tt statistics for the short and long timescales, respectively. This new distribution is well 
fitted with a power law of index ~ — 1 over the first 300 Myr (green). The purple dots 
are the data points from Fall, Chandar, & Whitmore (2005), as in Figure II. 5 (page 60). 
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pericenter passages: longer time intervals include the pre- and post-interaction history, 
when the tidal field is derived from an equilibrium configuration. 

When considering the possible destructive effect of the compressive mode, our new 
combination of the lus- and tt distributions can reflect both the formation and destruction 
rates of the clusters. That is, it is possible in a first approximation to identify the duration 
t of the compressive events with the age r of the clusters. In that case, the distribution 
derived in Figure III. 1 gives the number of clusters per age bin, and thus can directly 
be compared to the observational data. Surprisingly, the slope of our power law is very 
similar to the one obtained by Fall, Chandar, & Whitmore (2005) from observations of 
the Antennae, who found r^ 1 (recall Section II. 1.3, page 58). As a result, the duration of 
the compressive modes matches the age function of the star clusters, up to ~ 300 Myr. 
For older ages, our interpretation misses the secular evolution of the clusters due to 
internal phenomena, that rule their dissolution rate, in particular when the galaxies are 
quasi-isolated (i.e. for t > 300 Myr). 

Even if the link between our proxy and the empirical r _1 distribution is not a full 
chain of reasoning, the striking similarity of the two is still noteworthy. For the age 
interval considered here, taking into account the gravitational effect of the compressive 
tidal mode can account for the formation and dissolution rates of the clusters observed in 
the Antennae galaxies. Therefore, it is very likely that the evolution of the tides (going 
through a series of on/off compressive modes and shocks) influences the shape of the 
mass function and the dissolution rate of the clusters. As a consequence, the tides play 
a role in the empirical concept of "infant mortality" . 



1.3 Strength of the tides 

Up to now, we have mainly concentrated on the binary character of the tidal field 
(extensive or compressive), i.e. the very existence of the compressive modes. However, 
to compare its energy input with other mechanisms, one needs to quantify the strength 
of the tides for the various regions of the galaxies, at different steps of the merger. This 
will help to determine in which situation (where in the galaxy and for what duration), 
the tidal field is the most important, on the cluster scale. 

Figure III. 2 shows the distribution function of the the maximum eigenvalue of the tidal 
tensor from our Antennae simulation, at three different epochs. At any stage in the evo- 
lution of the merger, this distribution is by far not symmetric with respect to zero (which 
bounds fully compressive and extensive modes). Indeed, most of the particles have exten- 
sive (A max > 0) yet weak tides, while the smaller compressive mass fraction shows larger 
(negative) values. When the two progenitors are still well-separated (t = —100 Myr), the 
compressive part of the distribution (left-hand wing) contains a relatively small fraction 
of the mass. This fraction increases however during the interaction and reaches highly 
negative values. The merger clearly acts on the distribution by "expanding" its wings to 
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very negative values and, in a similar but milder way, to large positive values (extensive 
modes). Knowing that a given particle goes through a series of compressive/extensive 
switches (recall Figure 11.30, page 95), one may expect fast transitions from the left-hand 
side to the right-hand side of this distribution, and vice versa. During the major steps 
of the interaction, the dynamical time of the merger is shorter, which corresponds to a 
rapid motion of the clusters along their orbit (for those in the central region). That is, 
the switch from compressive to extensive (or the opposite) is more rapid and frequent 
than in the quiescent periods. The expansion and flattening of the distribution of the 
maximum eigenvalue however, suggest that the amplitude of theses modes can be higher 
at the times of the peri center passages. As a result, the tidal modes would be shorter 
but more intense than during a quiet phase. 

Figure III. 3 completes the picture by showing the distribution of the maximum eigen- 
value along the y-axis space coordinate. The Figure helps to understand that the ex- 
tremely compressive regions are found mainly in the nuclei, while the outer volumes yield 
compressive tides of amplitudes comparable to those of extensive modes. Note that only 
the maximum eigenvalue is shown here, which means that the minimum value associated 




Figure III. 2: Distribution function of the maximum eigenvalue of the tidal tensor of the 
Antennae model before interaction (blue), at the first passage (green) and during the 
merger phase (red). Negative values denote a compressive mode. 
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with a compressive mode (A min < A max < 0) has an even higher absolute value. 



1.4 Energy 

Knowing the strength of the field, we are now able to quantify the energy input by 
the tides. One can compare the role of the tidal field with that of other phenomena, via 
an estimate of the energy involved. Let's consider a cluster of mass M embedded in an 
isotropic tidal field, characterized by the unique eigenvalue A. The associated energy is 
given by: 

1 f M 

-Etides = --A / r 2 dm 
* Jo 

= - l -\aM R 2 t , (III. 7) 

where R t is the radius where the density of the cluster drops to zero, and a is a di- 
mensionless quantity enclosing the details of the mass distribution (see Section E.1.3, 




Figure III. 3: Left: spatial distribution of the particles in compressive (red) and extensive 
(blue) mode, in the orbital plane at t — 150 Myr. Right: absolute value of the maximum 
eigenvalue as a function of y. 
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page 234). The resolution of our simulations (M = 10 5 M and R t = 10 2 pc) corre- 
sponds well to the typical mass and size of a star cluster (Meylan & Heggie 1997). One 
can assume a certain mass distribution and set a (generally of the order of unity, see 
Section E.3, page 235) in Equation III. 7. As seen in Figure III. 2, the typical values for 
A are ~ 10~ 30 s~ 2 , which translates to a tidal energy of ~ i0 48-49 er g This does not 
compare well with the kinetic energy released by a type-II SN blast (10 51 erg) which 
would, therefore expand similarly with or without a tidal field. The winds of O-type 
stars may, however, be affected. Several authors evaluated the energy related to this 
process at the level of 10 46-48 erg (Cappa & Benaglia 1998; Martin, Cappa, & Benaglia 
2008), i.e. comparable or perhaps smaller than the energy of the tides. That is, a gas flow 
driven by massive stars would be slowed down significantly by the external effect of the 
tides and probably even kept bound within the cluster. In this case, the mass-loss due 
to stellar winds would be considerably reduced, which may prevent the early dissolution 
of the young cluster. 

In the TDGs region (recall Figure II.25.e, page 88), the compressive tidal mode acts 
on much shorter times but over an extended region. By applying the same estimates from 
Equation III. 7 to such a volume (M = 10 9 M and R t = 10 kpc), one gets 10 56-57 erg, 
i.e. much more than the energy release of a SN event. This value can be compared with 
a (rough) estimate of the number of type-II SNe that could form in such a region. Let's 
assume an initial gas mass of M gas = 10 9 M , converted into stars with a SFE of e ~ 10%. 
IMF | Then, consider a Salpeter (1955) initial mass function of the form £(m) oc m~ 2 35 , and 
get the mass fraction associated with the massive stars (m e [10, 120] M ) which will 
die as type-II SNe: 



Note that the slope of the Salpeter IMF for the high mass end is flatter than the one 
proposed by more recent works. For instance, Kroupa (2001) suggested £(m) oc m -2 7 
for M > 3 M , which would lead to less SNe and thus a smaller energy release into the 
ISM (see also Scalo 1998). That is, the Salpeter IMF gives an upper limit on the number 
of SN events, thus setting the maximum energy level to compete with the tidal field. 

With a mean mass of m = 22 M for the SNe, one has M gas e / m ,SNe-n / M « 10 6 
type-II SNe formed in the region considered. The associated instantaneous energy release 
is therefore 10 57 erg and compares well with the energy of the tidal field 5 . In other words, 
the compressive tidal field can compete with the usual stellar feedback (winds for the 
smallest star cluster forming regions, and supernovae explosions for the bigger TDG 
forming volumes) and alter the dynamical and chemical evolution of these regions, by 

b This estimate misses the details of the SFE and the IMF, as well as the actual effect of the feedback 
(not simultaneous in time and space for all SN blasts) on the ISM (turbulence, external pressure ...). 





(III.8) 



m £(m) dm 
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changing their binding energy. 



2 Star formation 

In this Section, we take one step further and include the tidal field in the chain of 
processes that form stars. To understand its role on the microphysics, one must recall 
the general concept of star formation. Stellar associations are formed by the collapse of 
a molecular cloud, when the contracting effects (gravitation, external pressure, shocks) 
overwhelm the internal pressure. In a first approximation, this can be encapsulated in 
the Jeans formulation summarized in Box III. 2. 



■ Box III.2 



The Jeans criterion describes the stability of an homogeneous sphere which 
is infinitesimally compressed from a density po to a new density p\ . By solving 
the fluid equations and taking the gravitational potential <f> into account, one 
gets the wave equation 




- v 2 V 2 pi - 47rGp pi = 0, 



(III.9) 



dt 2 



where v s is the speed of sound in the medium (see Binney & Tremaine 1987 for 
the detailed derivation). From this follows the dispersion relation between the 
frequency u and the wavenumber k of the solution pi(x, t) oc exp {i[kx — uit]): 



u 2 = v 2 k 2 - 47rGp, 



(111.10) 



That is, the growing solution for the perturbation exists if u 2 < 0, corre- 
sponding to an instability. Therefore, the compressed sphere is unstable if 




(111.11) 



or, in term of wavelength 




(111.12) 



Gp 



Aj is called the Jeans length. 
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Figure III. 4: Fragmentation of a molecular cloud (left) into proto-stellar cores (red, right) 
that accrete the surrounding material. The situation on the right-hand side snapshot is 
reached after two free-fall times. (Adaptation of Figure 1 of Klessen, Burkert, & Bate 
1998.) 



This corresponds to the Jeans mass 

Mj = ¥"°(y) 3 ' (IIL13) 

which defines the upper limit for objects stable to gravitational collapse. In 
other words, a gas cloud more massive than Mj has a proper gravitation 
that makes it collapse. In the first approximation of a homogeneous and 
non-turbulent medium, the Jeans criterion can be used to determined the 
conditions of the collapse of molecular clouds, i.e. the early phases of star 
formation. 



At smaller scale, the molecular cloud fragments into small clumps (see Figure III. 4) 
that continue to accrete the gas reservoir and cool through radiative transfer (Klessen, 
Burkert, & Bate 1998, Klessen & Burkert 2001). Turbulence and magnetic field play 
important roles in this process (see Shu, Adams, & Lizano 1987 and Mac Low & Klessen 
2004 for reviews), but are out of the scope of this document. 

In practice at cluster scale, the numerical simulations mimic the Schmidt-Kennicutt 
empirical law (recall Box II. 2, page 64) by setting a star formation rate depending on 
the gas (volume) density 

p* oc p g ' 5 if p g > p c and p* = else, (III. 14) 

p c being a critical gas density (see also Springel 2000, Springel & Hernquist 2003 and 
Bate, Bonnell, & Bromm 2003). Additional processes complete the model by reg- 
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ulating the star formation. In particular, stellar (SN) feedback is transferee! to the 
interstellar medium either via small kicks in the kinetic energy of the surrounding par- ISM 
tides in smoothed particle hydrodynamics (Navarro & White 1993), and/or through a SPH 
heating of the medium (Springel 2000; Kim, Wise, & Abel 2009). The cooling process 
can also be included, often thanks to the functions of Sutherland & Dopita (1993) for 
collisional ionization equilibrium. 



2.1 Galactic to cluster scales coupling 

During its collapse and fragmentation, a molecular cloud does not keep a regular 
shape but rather a filamentary structure including dense seeds (recall Figure III. 4). At a 
slightly more advanced stage, the young cluster will probably have a spheroidal profile, 
before getting a spherical shape, typical of globulars. This means that a young cluster 
is far from being spherically symmetric. Therefore, its host galaxy creates a torque that 
can modify the internal dynamics of the cluster. The goal of this Section is to describe 
this torque which couples the galactic scale to cluster-size elements. 

This problem strongly relies on angular dependencies for both the cluster and the 
galaxy. To ease the mathematical derivations in the general case, it is preferable to use 
spherical coordinates, as opposed to the more usual cartesian base. Hence, our idea 
is to decompose the quantities required to compute the torque in the base of spherical 
harmonics . In addition, Anderson (1992, see also Kawai & Makino 1998) provides a 
way to evaluate a potential at the position r from the value <f)(as) it takes on the surface 
of a sphere, called Anderson's sphere, of radius a: 

^ = ^ / E^+^G"^^)^ 8 )^ ( IIL15 ) 

for r < a, and 

= T / E( 2 ™ + l ) (-) p ™ (— ) ^ as ) ds ( IIL16 ) 

71 n=0 r r 

for r > a, where P n denotes the n-th Legendre polynomial. This allows us to compute 
the potential of either the galaxy or the cluster at any point, from its values on the 
Anderson's sphere, given by the spherical harmonics decomposition. (A step-by-step 
derivation is available in Appendix D, with conventions and notations.) 

For example, let's compute the torque created by a point-mass galaxy on a cluster. 
One needs the force of the galaxy on the stars of the cluster, i.e. minus the gradient of 

c Interestingly, the historical origin of spherical harmonics comes from Pierre-Simon de Laplace who, 
in Mecanique Celeste (1782), computed the gravitational potential of a set of point-masses, using the 
work of Adricn-Marie Legendre. The name "spherical harmonics" first appeared in Treatise on Natural 
Philosophy (1867) by William Thomson (aka Lord Kelvin) and Peter Guthrie Tait. 
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galaxy 



Figure III. 5: A point mass galaxy creates a torque on a star cluster. Setting the An- 
derson's sphere around the cluster eases the computation of either the potential of the 
galaxy at the positions of the stars, or the potential of the cluster at the position of the 
galaxy. 

the potential of the galaxy, at the positions of the stars. This can be done through the 
evaluation of the potential of the galaxy on the Anderson's sphere which includes the 
entire cluster (see Figure III. 5). 

Therefore, the torque created by the galaxy on the cluster is 



with p c the density of the cluster, and 0g the potential of the galaxy given by 
Equation III. 15, the stars being inside the sphere. 

Note that, because of Newton's law of reciprocal actions, it is also possible to compute 
the same torque from the potential of the cluster and the density of the galaxy (as in 
Appendix D): 



This time however, one needs to evaluate the potential at the position of the galaxy, i.e. 
outside the Anderson's sphere, and thus it must be computed thanks to Equation III. 16. 
In the case of a point mass galaxy, the integral reduces to the value of the integrand 
at the position of the galaxy. For more complex distributions, one can decompose the 
density pq in spherical harmonics on a shell centered on the cluster, but bigger than the 
Anderson's sphere. In fine, the problem is reduced to the computation of the torque 
created by one shell on another, concentrical one. The analytical solution in the general 
case is given in Section D.3 (page 229). 

To summarize, there are two ways of solving this problem: from the decomposition in 
spherical harmonics of 




(111.17) 




(111.18) 
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• the density of the galaxy and the potential of the cluster, 

• the density of the cluster and the potential of the galaxy. 



Configurations with a high degree of symmetry have been used as test cases of the 
formulas and their numerical implementation. The next step is to use the mass distri- 
butions obtained from our simulations of merging galaxies, and from simulations of star 
clusters. 

In the logical continuation of our numerical exploration, we choose to evaluate the 
potential of the galaxy on a sphere of similar size as the cubes used in the computation 
of the tidal tensor (recall Chapter I). 



Box III.3 

Technically, the potential of the entire galaxy is computed on a spherical 
grid thanks to getgravity. These values are passed to the tool S2Kit (Healy 
et al. 2003) which computes the coefficients of the spherical harmonics de- 
composition. Then, our code applies Equation D.60 (page 230) by summing 
on all the particles of the cluster to get the total torque created by the galaxy. 

Note that when the spherical grid is centered on a particle of the galaxy 
model, this central particle yields only a monopole term for the potential, and 
thus no torque on the cluster. Therefore, there is no need to subtract the 
value of its potential, by contrast with our method for the tidal tensor (recall 
Box 1.2, page 39). 

In practice, the sum on the spherical harmonics can be truncated to the 
sixth order, the coefficients of higher orders being negligible. 



To generate an anisotropic distribution of the mass of a proto-cluster (as in 
Figure III. 4), we set up configurations based on the cold-collapse of spherically sym- 
metric systems. A small number of particles (~ 10 3 ) is deliberately chosen to obtain 
a Poisson noise that favors the formation of substructures like clumps. Note that the 
turbulence modes present in fragmenting molecular clouds (see e.g. Mac Low & Klessen 
2004 and references therein) yield a more complex noise spectrum. 

First, a homogeneous sphere made of ~ 2000 particles is created and left with no 
internal velocity. Its cold-collapse rapidly gives out clumps of matter, as shown in the 
left-hand side of Figure III. 6. Second, a Plummer sphere is truncated to 90% of its 
total mass (~ 2000 particles too) and also undergoes a cold-collapse (right-hand side 
of Figure III. 6). Note that the density profile of the Plummer sphere implies a more 
centrally dominated collapse, the outskirts being more slowly affected than the core. 
Because of orbit crossing, the structure obtained is more elongated and resembles a bar. 
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These two transient structures have degrees of asymmetry close to those seen in the 
simulations of the collapse of a molecular cloud and thus, are expected to yield com- 
parable torque when embedded in a host galaxy. To compare their response to the 
scale coupling described above, Figure III. 7 plots the norm of the torque exerted by a 
point-mass galaxy situated at 10 kpc and 10 6 times more massive than the cluster. The 
corresponding statistics are given in Table III. 1 . Clearly, the bar mode (ex-Plummer) 
yields a stronger torque than the clumps (ex-homogeneous) and has a higher dispersion. 
This is explained by a weak torque when the bar is aligned with or perpendicular to the 
galaxy, and a much higher value for an oblique inclination, while the clumps always give 
out a comparable torque. 



o 
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Figure III. 6: A homogeneous distribution (top-left) forms clumps of matter (bottom- 
left) during its cold-collapse. The same transformation applied to a truncated Plummer 
sphere (top-right) yields a bar (bottom-right) . 
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Figure III. 7: Map of the torque created by a point-mass galaxy on the mass distributions 
shown in Figure III. 6 as a function of their orientation (8, ip) with respect to the galaxy. 



Table III. 1 : Statistics of the torques 



Model 




mean 


standard deviation 




[xlO 40 kg m 2 s~ 2 ] 


[xlO 40 kg m 2 s- 2 ] 


Clumps 


(ex- homogeneous) 


1.1 


0.4 


Bar (ex- 


Plummer) 


4.5 


1.8 
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Box III.4 

The effect of a torque can be compared with the auto-gravitation of the 
proto-cluster, in terms of energy. An order of magnitude estimate of the 
equivalent kinetic energy of the torque is given from the time derivative of the 
angular momentum: 

dL 
dt 

^m^, (111.19) 



At 

where (x) denotes the mean value of x, so that 



2 



-^torque ~ ^m(v) 



2{r)m 



(111.20) 



If one assumes that the galaxy creates a torque T = 10 40 kg m 2 s -2 during 
At = 1 Myr on a proto-cluster of mass m = 10 5 M and half-mass radius 
(r) = 5 pc, the equivalent energy is -Etorque ~ 10 44 erg. For comparison, the 
potential energy of the same proto-cluster would be 

Q « ?7T ~ 105 ° er §- ( IIL21 ) 
2(r) 

This means that a typical torque due to a galaxy on a proto-cluster does not 
compete with the gravitational energy and cannot modify its binding energy 
in a significant way. The orbital structure however, could be modified. 



To conclude, a point-mass galaxy can induce a stronger torque on an elongated proto- 
cluster than on a clumpy one. Such a torque may create or enhance the internal rotation 
of the cluster and thus reorganize its population. Note however that the torques com- 
puted are instantaneous and that the evolution of the proto-cluster, as well as the one 
of its environment, would change the orientation and/or the distribution of the matter, 
and thus the net effect of the torque. 

To go further, the point-mass galaxy is replaced by our model of the Antennae, pre- 
sented in Chapter II. Its potential is decomposed on an Anderson's sphere of radius 
a = 220 pc (i.e. the size of the cubes used to compute the tidal tensor) placed at the 
beginning of the compressive episode along the orbit of the particle D (next to the green 
x sign on Figure 11.30, page 95), i.e. 305 Myr after the first pericenter passage of the 
progenitors and in the western loop. This corresponds to a star forming region observed 
in the Antennae galaxies. For the cluster, the distributions obtained by cold-collapse are 
scaled to get a total mass of 2 x 10 5 M and a mean radius of ~ 5 pc. 
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Figure III. 8: Map of the ratio of the torque created by the Antennae and a point-mass 
galaxy on the mass distributions shown in Figure III. 6. See text for details. 



The main source of gravitation is now much closer and more massive than in our 
point-mass example: ~ 1 kpc instead of 10 kpc and twice the mass (corresponding to 
two nuclei). As a consequence, the torque computed is stronger, in average, by a factor 
~ 200. To have a fair comparison of the values of the torques, the axes of reference 
for 9 and <p have been chosen so that the dipole of the potential of the Antennae (i.e. 
the direction of minimum potential) has the same angular position than the point-mass 
galaxy in the previous example. 

Because the torque partially depends on the mass distribution of the cluster (e.g. 
possible alignment of the bar with the galaxy), we find similarities in the angular map of 
the torque created by the Antennae and in Figure III. 7. The ratio of the two is plotted 
in Figure III. 8. As estimated before, the mean values of the ratio are 207 and 203 for 
the clumps and the bar respectively, and the standard deviations reach 140 and 128. 
This shows that the torque on the bar depends slightly less on the mass distribution of 
the galaxy than the one on the clumps. This can also be interpreted as a more robust 
response of an elongated structure in a varying potential (orientation and strength), e.g. 
along the orbit of a cluster-size element. 

It is possible to estimate the equivalent energy of these torques. The proto-cluster 
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being much closer to the galaxy than in the point-mass experiment, the time during which 
the torque may be effective is much shorter. Using an estimate from the circular velocity, 
one finally gets an energy ~ 5 x 10 2 times higher than the one obtained for the point 
mass. However, these estimates have a strong dependance with, e.g. the distance and 
mass of the galaxy. Therefore, such orders of magnitude remain imprecise. For a finer 
description of the coupling between the galactic and cluster scales, numerical simulations 
are required. In this matter, a star cluster must be simulated using the external effect of 
the galaxy, as we do in the Section 3. 



2.2 Star formation efficiency 

The existence of a tidal field introduces a new term in the expression of the virial 
theorem. (Appendix E details its derivation.) Therefore, the criterion for the survival 
of a young cluster after its gas expulsion is changed. To illustrate this, we repeat the 
analytical work of Hills (1980) and Boily & Kroupa (2003a) made for clusters in isolation, 
in the case of an external tidal field (see Section E.4, page 238). 

In isolation, a young cluster dissolves if the velocity dispersion of its stars is not 
balanced by the gravitational well. If it is the case, the cluster is not able to keep its 
material bound. When the very first stars form in the cluster, their velocities are set 
according to the gravitational potential, built by both the stars and the gas. When the 
gas spills over to the galactic ISM, the stars remain with the same velocity dispersion, 
while the potential well has become shallower. In other words, the escape velocity is 
reduced to a level possibly below that of individual stars. This depends on the fraction 
of the mass lost, i.e. the mass fraction of the gas, with respect to the total mass of the 
cluster. This quantity comes from the SFE e: 



Hills found that the cluster survives if e > 50%, i.e. if at least half of the initial mass is 
converted into stars (see also Section E.4.1, page 238). 

When the same proto-cluster is embedded in a tidal field, the dynamics is changed. 
Again in Section E.4 (page 238), we have extended the previous derivation by taking 
into account the additional effect of an isotropic tidal field (i.e. represented by an unique 
eigenvalue A), starting from the work of Fleck (2007), and extending his results to other 
transformations. 

In this case, the response of the cluster to the tidal field is a strong function of its 
density profile and thus, the calculation is more involved. It is very likely that this 
distribution changes after the gas expulsion. Unfortunately, a general solution does not 
exist. For simplicity, let's consider a proto-cluster which will undergo a homologous 



Mi 



gas — 



M-M* 
M(l-e). 
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transformation, i.e. keeping the same shape but not the same amplitude d of density p Q . 
We note A the ratio of tidal to gravitational energies: 

Xa R 2 r 

A = ^ (....23) 

(see Appendix E for the definitions and notations). Equation E.46 (page 241) tells us 
that, after an instantaneous expulsion of the gas, the cluster remains stable if 



I | M - 1.4 

rv 



a' (R\ 

1 + - IT 
a V Rt 



< 0. (111.24) 



By homology, we have a = a' and r' v /r v = R' t /R tl which leads to 

e R't ( 1 „ A „ . x '" 



2 + afU" e_M >> +M W - a (IIL25) 

This relation constrains the radius i?' t of the cluster after the gas expulsion, once one 
knows (i) the SFE e and (ii) the strength of the tidal field, embedded in the term A oc A. 
Figure III. 9 gives some examples of expansion of the cluster, for several tidal fields. In 
extensive mode, a solution (i.e. a positive radius R't) exists only for a certain range of 
SFEs. Note that, without a tidal field (A = 0), one retrieves the solution of Hills (1980) 
and the e > 50% limit. However, the most important conclusion is that the presence of 
a compressive field allows the cluster to reach a new equilibrium stage after the gas loss, 
for all values of e. 

In other words, when the gas is expelled, the kinetic energy of the stars is too high with 
respect to the rest of the energy to keep the cluster bound: it is supervirialized. Therefore, 
it must expand. Our analytical solution gives the range of possible expansions (y-axis on 
Figure III. 9) that keep the cluster bound. In an extensive mode, a too severe mass-loss 
(low e) leads to the dissolution, while a higher SFE allows some "bound" expansions. 
However, if the tidal field is strong (large A), important expansions correspond to the 
dissolution of the cluster due to the extensive mode. By contrast, in compressive mode 
the tidal field cannot dissolve the cluster. Therefore, the expansion has only a lower 
limit, set by the kinetic energy of the stars. 

As seen in Figure 11.30 (page 95), a young cluster will likely experience a switch from 
compressive to extensive mode. If it forms in a compressive, supervirialized cocoon and 
then goes to extensive mode, it will likely dissolve. (This is visible on Figure III. 9 where 
the compressive mode offers more possible configurations than its extensive equivalent.) 
Such a process plays an important role in the "infant mortality" process, and thus affects 
the death rate of the clusters and their age distribution. Note that it is not straightfor- 
ward to conclude on the mass function, because the link between the tidal energy and 
M is only known once the mass profile of the cluster is set. 



d Keep in mind that is assumption is very rough: the external layers of the cluster are much more 
affected than its core which inevitably leads to a reorganization of the matter, and thus a change of the 
density profile. 
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2.3 Multiple main sequences 

In addition to possible mass loss, the repeated on/off evolution of the compressive tidal 
mode can influence the chemistry of the cluster. In this Section, we introduce an idea 
that may (partially) explain the existence of observed features in the color-magnitude 
diagram of several massive clusters. The reader should keep in mind that the following 
arguments are purely qualitative and speculative. The exploration of this point requires 
a detailed numerical treatment that deserves a study of its own. 

In our common understanding of star clusters, we assume an initially homogeneous 
chemical composition and a coeval formation of the stars. This idea is often referred to 
as the SSP. However, recent observations of resolved massive star clusters have revealed 
multiple populations (see e.g. Bedin et al. 2004; Piotto 2009; Milone et al. 2009b) 
in several cases, denoting a more complex star formation history and/or chemistry at 
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Figure III. 9: Expansion of a cluster after the expulsion of the gas, in an extensive tidal 
field (A > 0, blue curves), without tides (A = 0, green) and in compressive mode (A < 0, 
red), as a function of the star formation efficiency e. The numerical value on each curve 
is A, i.e. the ratio between the tidal and gravitational energies. The areas filled with 
red and blue dots correspond to bound configurations for clusters, for A = —0.01 and 
A = 0.01, respectively. 
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the time of formation. These appear in color-magnitude diagrams as multimodality in 
sequences or branches, as visible in Figure III. 10. 

The first discovery of such a multiple main sequence has been made in u Cen (the | MS 
most massive cluster in the Milky Way), which displays a split into a dominant red and 
a secondary blue MSs, at all radii. Spectroscopy shows that the bluer sequence has twice 
the light elements abundances (He, CNO) of the red one (Piotto et al. 2005). These 
peculiarities confirmed that u Cen is an anomalous cluster, possibly the remnant of a 
dwarf galaxy (see Lee et al. 1999). However, the presence of multiple populations of stars 
has been detected in other, perhaps less peculiar, clusters in the Galaxy like NGC 2808 
(Piotto et al. 2007), NGC 1851 (Milone et al. 2008) or others (not yet published, see 




Figure III. 10: Color-magnitude diagram of NGC 1851 made with HST ACS by Milone 
et al. (2008). Two SGBs are clearly visible in the upper part of the diagram. (Data 
points from Milone et al. 2008.) 
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Piotto 2009 for a list of preliminary results). Recently, clusters showing similar features 
in their color-magnitude diagrams have been found in the Small and Large Magellanic 
Clouds (Glatt et al. 2008), where it has even been mentioned that most of the clusters 
SGB | contain multiple populations (detected via splits in the turn-off and sub-giant branch of 
their color-magnitude diagrams, see Milone et al. 2009a and Figure III. 10). 

The explanation of such observations is still debated. Because each cluster yields 
different properties of its branches, it is likely that the origin of the multiple populations 
is far from being unique. A first possibility would be the formation of new stars during 
the encounter between the cluster and a giant molecular cloud. Bekki & Mackey (2009) 
have shown that this scenario would explain the presence of two turn-offs, one due to the 
initial cluster and the other from the stars formed during the interaction with the cloud, 
at a different epoch. 

Alternatively, Pflamm-Altenburg & Kroupa (2009) suggested that a massive cluster 
could accrete some of the gas met along its orbit through the galaxy. Again, this new 
gas would be processed into a new generation of stars, explaining both the age and 
abundances differences. 

Another possibility involves the merging of two distinct star clusters. One can imagine 
that, in binary systems of clusters, a loss of angular momentum could lead to a merging of 
the two populations of stars into a single structure. However, the age difference between 
the two remains difficult to explain because both clusters are supposed to share the same 
origin, and thus to have similar properties, like the age of their stellar content (see e.g. 
Dieball, Miiller, & Grebel 2002). 

These ideas are based on the formation of stars from an external source of gas. By 
contrast, one can imagine the creation of a secondary population from an internal source 
of gas, i.e. the miracluster medium itself. This implies to be able to retain the ejecta 
of the stellar feedback (winds and SN blasts) within the cluster and introduces a "self- 
pollution" . The enrichment of the ISM within the cluster by the first stars would be the 
source of a new generation, with higher CNO (carbon, nitrogen, oxygen) abundances, 
as observed. The main issue of such a scenario is that the ejecta have to remain in the 
cluster instead of being blown away. Several theories suggest to consider the matter from 
the winds of asymptotic giant branch stars (D'Ercole et al. 2008) or massive rotating 
stars (Decressin et al. 2007; Decressin, Baumgardt, & Kroupa 2008) whose winds are 
slow. After the violent gas removal due to type-II SNe, these ejecta could gather in the 
cluster center and be transformed into stars. 

In this manner, the presence of a compressive tidal mode could be the missing link that 
helps to bring and/or retain the matter within the cluster. Indeed, if the cluster moves 
through an external source of gas (a molecular cloud or a second cluster), the compressive 
field could increase the binding energy of the main cluster so that this external amount 
of matter is captured and retained. Similarly, the ejecta from the cluster members would 
be more easily kept bound. 
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On the one hand, the accretion from an external source would not be specific to certain 
chemical elements and thus, even heavy ones (e.g. iron) should be detected in the second 
generation of stars. With the retention of ejecta on the other hand, the new abundances 
should be more similar to those of the first, original generation. 

The fine details of these processes (winds, accretion, enrichment, impact of the tides) 
are not well understood yet and it is too early to pin down the exact effect of the tidal 
field on such a complex star formation history. 



2.4 Schmidt-Kennicutt and the compressive mode 

As mentioned in Box II. 2 (page 64), the SFR can be empirically linked to the power 
w 1.4 of the surface density of the gas, via the Schmidt-Kennicutt relation. In addition, 
our results suggest to connect the formation of stars to the concept of compressive tidal 
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Figure III. 11: Surface density of the particles in compressive mode S c compared to those 
of all the particles S. The red dashed line marks the equality of the two quantities. 
The blue line shows the best fit of the upper envelope of the scattered distribution, as a 
power-law of index w 1.40. The apparent binning is due to an integer number of particles 
per pixel. 
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mode. Therefore, it seems natural to seek a relation between the surface density of the 
gas and the compressive mode. To explore this further, we plot in Figure III. 11 the two 
surface densities £ (all particles) and £ c (compressive particles only) for all the pixels of 
the map presented in Figure 11.27 (page 90). 

Because one cannot find more particles in compressive mode than the total number of 
particles, the line £ c = £ is the logical upper limit. We define the envelope of the points 
displayed here as the maximum value of E c taken in bins of 20 M Q /pc 2 . It can be well 
adjusted with a power-law E c ex E L4 ° (with the standard deviation a = 26.2 M /pc 2 ), 
i.e. the same expression as the empirical Schmidt-Kennicutt relation. This means that 
the surface density of compressive modes is limited by the same relation that links the 
SFR with the surface density of the gas. 

Note however that our relation E c oc E L4 ° concerns the envelope and thus is a limit: 
many points stand below this line on Figure III. 11. Such a dispersion does not exist 
in the Schmidt-Kennicutt relation, where all the data points are close to the empirical 
law. All together, these results suggest that (i) again, the compressive mode plays a 
role in the formation of star as it shows a similar relation with the density as the SFR 
does and (ii) the large dispersion in our relation appeals for a finer prescription for 
star formation, including sub-scale physics and hydrodynamics. Such a more detailed 
comparison between the evolution of the compressive modes and the SFR throughout 
the interaction history is feasible when considering hydrodynamical runs. 

2.5 Hydrodynamical studies 

In terms of mass fraction as well as intensity, the compressive tidal field is clearly 
associated to the major events in the interaction history of the merger as revealed by 
Figure 11.23 (page 86). One can also identify the peaks around the pericenter passages 
with those noted by several authors who studied the SFRs from simulations of merging 
or interacting galaxies. For instance, di Matteo et al. (2007) presented a statistical 
study of 240 SPH simulations of interacting galaxies and explored the evolution of their 
SFRs (see also di Matteo et al. 2008). They noted that the pericenter passages could be 
identified with sharp increases in the SFR, by a factor depending on the fine details of 
the morphological types of the progenitors and the orbital parameters 6 . 

In our own description, the impact of tidal stripping during close passages is dimin- 
ished because a pressurized component is missing. In addition, we do not take into 
account the formation of stars at the first passage, that inevitably depletes the gas reser- 
voir, leaving a lower gas mass available to form stars during the second passage. However, 
even without the response of the gas to the interaction, the evolution of the mass fraction 
in compressive mode (i.e. the relative importance of this mode) matches well the curve 

c Note however that they suggested that not all interactions are to be systematically linked with a 
starburst activity. 
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of the SFR versus time. 

We have recently extended our palette of tools to explore the history of star formation, 
with two techniques that include hydrodynamics: adaptive mesh refinement and SPH. | AM R 
We have applied these tools to models of the Antennae, and we begin to compare their 
results with the properties of the compressive tidal mode. 



With AMR 



An hydrodynamical study of our model of the Antennae is currently conducted in 
collaboration with researchers at the CEA f (PI: D. Chapon). With the initial condi- 
tions presented in Table III (page 69), a gaseous disk of 10% of the baryonic mass has 
been added to both galaxies. Gravitational and hydrodynamical equations have been 
solved thanks to the AMR code Ramses (Teyssier 2002) with 14 levels of refinement 
(Figure III. 12), which gives a maximal resolution of 12 pc. Figure III. 13 shows the mor- 
phology of the system at t ~ 100 Myr, i.e. between the two pericenter passages. Because 
of the additional gaseous mass and its response to the interaction, a small difference 
exists in the timing between the two passages: the AMR run gives out a separation of 
260 Myr instead of 300 Myr for the purely gravitational model. 



Box III.5 

The AMR technique allows to place the computational power where needed. 
The code starts with a regular box (in our case, a cube of 200 kpc) and 
checks for a refinement criterion (here, a certain density level). If the criterion 
is fulfilled, the cube is refined into cells, which are also examined for the 
same criterion and so on (in a similar way as a tree-code, recall Section 1.1.3, 
page 29) until either the criterion is not fulfilled, or the maximum level of 
refinement is reached. In this way, it is possible to have an extremely high 
resolution in regions of interest, like cluster forming zones, and to keep a low 
resolution in e.g. empty spaces to keep the computation fast. Figure III. 12 
illustrates this point. 

Once the mesh has been built, the equations of (magneto-)hydrodynamics 
are solved in each cell. In addition, the stellar particles are evolved thanks to 
a traditional Poisson solver for the total gravity (gas + stars + DM). 



In the preliminary runs, no feedback nor magnetic field have been introduced. The 
star formation recipe is based on a Jeans mass taken equal to 2 x 10 5 M (matching the 



Commissariat a l'Energie Atomique, Saclay. 
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Figure III. 12: Grid used in the simulation of the Antennae galaxies, viewed in the orbital 
plane at t ~ 100 Myr (i.e. as in Figure Il.lO.d, page 70). Black, green, blue and red 
correspond to the refinement levels 5, 6, 7 and 8. In total, the mesh counts ~ 10 6 cells (up 
to level 14). For graphical purposes, only the densest regions (in projection) are plotted, 
which explains that some cells seem to be not correctly split into four sub-elements in 
this Figure. 



124 



2 - Star formation 




Figure III. 13: Hydrodynamical simulation of the Antennae with the Ramses code. This 
snapshot corresponds to t « 100 Myr (i.e. as in Figure II.25.d, page 88). The density of 
gas is shown in green, red represents the old stars while blue denotes the newly formed 
stars. Thanks to the AMR, the largest pixels cover ~ 1.5 kpc while the smallest ones 
(e.g. the blue seeds of star formation in the bridges) represent 12 pc only. (Image from 
D. Chapon, private communication.) 
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t [Myr] 

Figure III. 14: SFR from an AMR simulation of the Antennae. See the text for details. 

typical value of the massive clusters) via a polytropic equation of state ff P oc p 7 . The 
formation rate of clusters is fixed at a high value, close to the one set by the numerical 
resolution limit, which gives a SFR that might be over-estimated. A deeper exploration of 
the parameters of the equation of state, the formation efficiency and the cooling function 
will be done in the next months. 

Figure III. 14 shows the evolution of the SFR during the simulation. Despite the 
orbital difference between the AMR and the dry runs, the SFR matches surprisingly well 
the trend seen in the mass fraction in compressive mode (thumbnail in Figure III. 15, 
below). Indeed, the first pericenter is clearly marked by a sharp peak, followed by a slow 
decrease up to the second passage. We note a difference however: the second major peak 
in the SFR is lower (by 20%) than the first one. In comparison, the two peaks are of 
similar amplitude for the mass fraction in compressive mode (dry run). This difference 
can be explained when considering the dynamics of the gaseous component. During the 
first interaction, the gas reservoir has been depleted by (i) tidal stripping and (ii) its 
conversion into stars. In other words, a fraction of the gas has been spread out (recall 

g When collapsing, a molecular cloud reaches a density so high that its heat cannot be radiated away 
(Mac Low & Klessen 2004). That is, it becomes opaque to radiative cooling and should no longer be 
considered isothermal (7 = 1). That is why the equation of state adopts the polytropic form P oc p 7 , 
with 7 > 1. 
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the Hi tails detected by Hibbard et al. 2001), and the rest has been processed to form 
the stars, with a certain efficiency e. These two reasons combined explain that a smaller 
amount of gas is available at the second passage to form new stars. That is why the 
second peak in the SFR is lower than the first one. For the case of the dry run on the 
other hand, the only difference between the two events is the mass-loss due to the creation 
of the tails (~ 30% of the total mass of the stellar disk). This should lead to a smaller 
second peak: however, because of the loss of angular momentum of the progenitors, the 
pericenter distance is much smaller (recall Figure 11.15, page 76), which enhances the 
fraction of compressive modes*. In turn, the combination of these two effects (mass-loss 
and close passage) balance each other so that the second peak turns out to be of similar 
amplitude to the first one. 

Even if our hydrodynamical study with Ramses is still at a preliminary phase and some 
tuning of e.g. the star formation recipe is required, we note a strong correlation with the 
dry run presented in the previous Chapter, and in particular between the evolution of 
the mass fraction in compressive mode and the one of the SFR. In addition to the good 
matches between the location of the compressive tidal field, their duration, their energy 
and the properties of the young clusters observed in the Antennae, this new relation 
with the SFR strongly suggests that the compressive mode is linked to the formation of 
substructures. An important work with our Ramses simulations will be done during the 
following months to fully explore and quantify the very first points we have highlighted 
here. One of the key question to be addressed pertains to the existence of a time delay 
between the creation of a compressive mode and the actual formation of proto-stellar 
seeds from the fragmentation of the ISM. 

Note that the AMR technique has already been used to study star formation in inter- 
acting galaxies by Kim, Wise, & Abel (2009). They set up a planar encounter of two disk 
galaxies and derived the evolution in time of the SFR, with a resolution of 2 x 10 3 M Q . 
There again, peaks exist at the times of the close passage, denoting a shock-triggered 
formation, as mentioned by Barnes (2004) with SPH runs. As a complement to this 
work, our new Antennae model with Ramses will provide a good sample to explore the 
details of the connections with the tidal field. 



With SPH 

In parallel with our AMR study, a brand new SPH model of the Antennae galaxies 
has been presented in Karl et al. (2010), using the code Gadget2 (Springel 2005). This 
model is part of a parameter study aimed at finding better initial conditions for the 
modeling of the Antennae galaxies. Note that the initial conditions differ from the purely 
gravitational run presented in Chapter II. In particular, the SPH model yields a longer 
orbital period. The resolutions however, are comparable. As visible in Figure III. 15, the 

h This point will be developed in details in the next Chapter. 
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SFR is, again, enhanced at the times of the two pericenter passages. We note that the 
first peak is much less pronounced than the second one, and that the associated "burst" 
of star formation is diluted in time over ~ 500 Myr. 

The reason for this is still under investigation but seems to be linked with the low spin 
of the disks which makes them more compact and thus, more resistant to the effects of 
the interaction. The gas is then slowly stripped and can fall back onto the disks during 
a much longer period, which translates into a diluted "burst" of star formation. In this 
case, the mild decrease of the SFR at ~ 400 Myr would correspond to the depletion of 
the intergalactic gas reservoir. For the second passage, the progenitors have already been 
affected which could mean a much more efficient and rapid formation of stars, as noted 
on Figure III. 15. 

Furthermore, the SPH model uses massive Hernquist DM halos while the dry run and 
the AMR model have lighter, isothermal halos. This additional difference may also play 
an important role in the formation of stars. The influence of the shape of the DM profile 



i — q 




-200 200 400 600 

t [Myr] 

Figure III. 15: SFR from the SPH simulation of the Antennae presented in Karl et al. 
(2010). Arrows indicate the two pericenter passages. The evolution of the SFR for 
t > 600 Myr is still being computed. The thumbnail shows the shape of the evolu- 
tion of the mass fraction in compressive mode from the dry simulation (as presented in 
Figure 11.23, page 86), reminder. 
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is studied in Chapter V in more detail. 

All the concepts and results presented so far show a correlation between compressive 
tidal field and star formation, in our models of the Antennae galaxies. It is important 
however to confirm that our conclusions do not exclusively apply to this merger but can 
also be extended to a broad range of configurations. Therefore, an exploration of the 
influence of morphological and orbital parameters on the characteristic of the compressive 
tidal modes is presented in the next Chapter, to allow a comparison with the trends in 
the SFR noted by e.g. di Matteo et al. (2008). 



3 Stellar dynamics 

Before investitaging other galaxies, let's focus on the early life of a cluster, i.e. after 
the formation phase. In addition to triggering or strengthening star formation, the host 
galaxy can modify the morphology and internal evolution of a cluster through the tidal 
field. Thanks to the tidal tensor, it is possible to take into account the potential of the 
galaxy on the spatially extended cluster in a relatively simple, yet realistic way (i.e. for 
non constant and non isotropic tides). In particular, the tidal field seen by a cluster along 
a non regular orbit (like those found in mergers) changes much more rapidly than the ones 
typically described in the literature (see e.g. Baumgardt & Makino 2003). Therefore, 
our method allows to explore the role of the large scales (galaxies) on the smaller scales 
(star clusters) for all kind of orbits. 

However, the study of the evolution of a cluster requires a more precise description 
than the collisionless iV-body approach we have followed so far. As presented in Box 1.1 
(page 28), globular clusters are collisional systems. Therefore, a numerical treatment 
of close encounters between stars, and binaries must be used. Furthermore, one should 
account for stellar initial mass function, stellar evolution, feedback... The code used 
in our galactic exploration (i.e. gyrf alcON) does not consider these aspects and thus, 
a different tool is needed. Nbody6 (see Aarseth 1999 and references therein) and its 
version for parallel computers Nbody6++ (Spurzem, Baumgardt & Ibold 2003) fit our 
requirements as they allow to simulate star clusters and the evolution of their stellar 
contents. 

Based on an idea of Fleck (2007), we have adapted Nbody6++ so that it can take 
an external tidal field into account via the associated tidal tensor. In this new code, 
called Titens, the initial tensor is strictly equal to zero and then, it is rapidly grown to 
reach the values required 1 . This ensures continuity and coherence of the initial setting of 
the cluster, which is built without external field. This growth is done over a couple of 
timesteps, i.e. ~ 10 4-5 yr. 

1 An immediate growth of the tensor would lead to a dramatic change of the energy of the system, 
which would make Titens detect an anomaly and stop the computation. 
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Box III.6 

Tit ens reads the nine values of the tidal tensor from the galactic simu- 
lation, plus the associated time, in an external file. The timestep of Titens 
(~ 10 4 yr) is much shorter than the one of gyrf alcON (~ 10 6 yr) and thus, the 
tensor is badly sampled for a direct use in a stellar dynamics run. Therefore, 
its values are interpolated when needed by Titens, using a linear scheme for 
the first snapshot and a quadratic interpolation for all the others. 

Note that both the time and the tensor must be scaled properly, because 
gyrf alcON and Titens do not share the same system of units. 



It is now possible to model a cluster with a given mass profile (e.g. Plummer, King, 
...), an IMF and stellar evolution (Hurley, Pols, & Tout 2000). For the preliminary runs 
presented here, we focus on the effect of the tidal field and thus, temporarily switch-off 
stellar evolution. Nonce, all the stars in our clusters have the same mass. Therefore, the 
effects of a mass function such as a non-trivial luminosity function and mass segregation 
do not occur. That is, the role of the tidal field only is made easier to detect and to 
analyze. Obviously, the full possibilities of Nbody6++ will be exploited right after these 
first studies. 

3.1 In an isotropic tidal field 

For the first examples, five isotropic and constant tidal fields have been artificially set, 
with values typical of what we found in the galactic runs: 

• an isolated case: A = 0, used ClS db reference, 

• a weak field (extensive or compressive): A = ±2 x 10 -30 s~ 2 , 

• and a strong field : A = ±6 x 10~ 30 s~ 2 . 



In addition, we have created a "transition" tidal field which is compressive for the first 
25 Myr and then becomes extensive forever, with the values of the weak regime in both 
cases. 

Note that an isotropic extensive tidal tensor corresponds to a negative density (its 
trace being positive, see Equation 12, page 11). Although not physical, such cases allow 
a simple analysis and a direct comparisons with ID simulations (e.g. Fokker- Planck, 
Monte-Carlo). The possible mass- loss due to the field will be, in this case, over-estimated. 

The effect of the tidal field depends on the density of the cluster. Therefore it is 
expected to be stronger in the outer parts of a mildly concentrated object than close to 
the center of a dense profile. To confirm this idea, King models with various parameters 
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(^o/ " 2 ) have been set-up: for a given total mass, a small parameter corresponds to a 
peaked profile and thus, a small half-mass radius, while a higher number is associated 
to a more extended cluster (see Binney & Tremaine 1987). Figure III. 16 sums up this 
by showing the surface density profiles of the three models we have used: King5, King7 
and King9. Every cluster is made of 10 4 stars of exactly 1 M each. The profiles do not 
change significantly over 100 Myr. The stability of the models in isolation is key in this 
study: it ensures that intrinsic phenomena (e.g. relaxation) do not bias the effect of the 
external field. 

The evolution of some Lagrange radii is shown on Figure III. 17 for the King7 model. 
The inner regions are not strongly affected by the tidal fields used, and the associated 
Lagrange radii (e.g. 50%) remain fairly independent of the strength of the tides. In the 
outer parts (70%, 90%) however, the extensive mode spreads efficiently the matter so 
that the "size" of the cluster increases. The compressive fields have a weaker, opposite 
effect on the radii but are still clearly distinguishable from the isolated case. 




Figure III. 16: Surface density profiles of the three King models in isolation, at initial 
time (solid curves) and after 100 Myr (dotted curves). Only a small evolution due to 
relaxation is visible. Vertical ticks in the upper part mark the half-mass radii of the 
models. 
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When embedded in the "transition" field (black curve), the cluster first gets a high 
potential energy: it becomes supervirialized. When the tides switch to extensive mode 
(t > 25 Myr), the potential gets flatter and the velocity dispersion of the stars is too 
high to maintain equilibrium. Therefore the outer regions are violently expelled from 
the cluster: this translates into the growth of the 70% and 90% Lagrange radii, much 
quicker than for the constant extensive field of the same strength (solid blue curve). 

On the one hand, the same effect is even more dramatic for the King9 model, as shown 
in Figure III. 18. In this case, the inner parts of the cluster are also affected and only 
the core (~ 10% radii) remains fairly stable. On the other hand, the King5 model is too 
dense to be really modified by the tides. Even the outer regions keep the same size than 
in isolation: no effect of the field is visible on the 90% Lagrange radius. 




I i i i i i i i i i i i i i i i i i i i I 

20 40 60 80 100 

time [Myr] 



Figure III. 17: Evolution of the 50%, 70% and 90% Lagrange radii of the King7 model in 
isolation (green), in a weak field (solid curves) and in a strong field (dashed curves). Blue 
corresponds to the extensive fields while red denotes the compressive modes. The solid 
black lines are associated to the "transition" field (i.e. weakly compressive for 25 Myr 
and weakly extensive afterwards.). See text for details. 
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Figure III. 18: Same as Figure III. 17, but for the 10%, 20% and 30% Lagrange radii of 
the King9 model. 



3.2 In the Antennae 

These examples help to understand the behavior of a cluster when embedded in a 
simple tidal field. However, the previous Sections and Chapters revealed that merging 
galaxies yield tides that are far from being constant or isotropic. Therefore, we repeat the 
same experiment, with the same clusters, but this time in tidal fields obtained from the 
simulation of the Antennae (see Chapter II). The orbit of the particle B has been chosen 
(recall Figure 11.30, page 95) because it has a long compressive episode induced by the 
first pericenter passage, and a moderated (yet significant) strength. The particle D gives 
out a comparable tidal field at the second passage, but with a stronger extensive regime. 
Figure III. 19 displays the evolution of the maximum eigenvalue of the tidal tensors of 
these particles, and the time lapses considered for the Titens runs. From now on, the 
time reference (t = 0) will be taken at the start of these lapses. 

Figure III. 20 shows the evolution of the Lagrange radii in the inner parts (30%, 50%) 
and in the outer regions (90%) along the two orbits, for the King7 and King9 models. 
As before, the King9 cluster is more deeply affected by the tidal field. The effect of the 
stronger extensive mode of the orbit D (with respect to those of the orbit B) is marked by 
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a more pronounced expansion of the external layers of the clusters. However the general 
behaviors remain fairly comparable along both orbits. 

The violent expansion of the outer parts of the King9 model implies a change of 
its morphology and thus, a modification of its visual appearance. To illustrate this, 
we now take an observational perspective. The surface density profiles of the clusters 
are shown in Figure III. 21. Because all the stars in our runs have the same mass, this 
directly corresponds to a surface brightness profile. However, at these positions inside 
the Antennae galaxies, the surface brightness surrounding the cluster is far from being 
null and should be taken into account. To do so, the average density at the position of 
the cluster (i.e. along the orbits of B and D) is computed, and an equivalent brightness 
for the background is determined. (The equivalent surface density is shown as dotted 
lines in Figure III. 21.) We set the detection level at twice the value of the brightness of 
the background. Thus, all the outer, diffuse regions of the cluster are not detected. 
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Figure III. 19: Evolution of the maximum eigenvalue of the tidal tensor of the particles B 
and D (see Figure 11.30, page 95). The x-axis indicates the time of the galactic simulation. 
The colored parts correspond to the fields extracted and given to Titens. Both start 
with a compressive mode (red) lasting 20 Myr, before going to an extensive regime 
(blue). In this last step, the field is stronger for the particle D than for B. 



134 



3 - Stellar dynamics 



Box III.7 

The brightness of the background depends on the column density along 
the line of sight, and thus, will change with the orientation of the galactic 
model. To avoid this, we compute the local volume density at the position 
of the cluster in the Antennae, thanks to the trace of the tidal tensor (recall 
Equation 12, page 11). Then, the volume density is integrated over the depth 
of the cube used to compute the tensor, and the surface density obtained is 
converted into a brightness, using the same scale as for the cluster itself. 



In the upper part of Figure III. 21, vertical ticks mark the half-mass and half-light 
radii: the half-mass radius takes into account all the stars of the initial cluster, while the 
half-light radius only considers the stars which have been detected. 

Figure III. 22 displays the surface brightness of the three King models embedded in the 
the tidal fields of the orbit B at t = 20 Myr (i.e. at the end of the compressive regime) 
and t = 40 Myr (in extensive mode). Again, the peaked King5 is almost not affected by 
the external field. Contrarily, the King7 and King9 models are strongly disturbed and 
their two radii increase. If one considers that a star which is not detected is in the field 




20 40 60 

time [Myr] 



Figure III. 20: Evolution of the Lagrange radii of the King7 (top) and King9 (bottom) in 
isolation (green), along the orbits of the particles B (solid red/blue lines) and D (dashed 
lines). The parts of the curves in red denote a fully compressive mode (first 20 Myr). 



135 



Chap. Ill 



Tides and star formation 



and not anymore in the cluster, the mass-losses at t = 40 Myr are ~ 8%, 29% and 76% 
of the initial mass (10 4 M ), for the King5, King7 and King9 models respectively. 

Along the orbit of the particle D, the models give out a comparable response to the 
tidal field (Figure III. 23). As already visible in Figure III. 20, the half- mass radii remain 
similar to those found in along the orbit B. However, the half-light radius of the King9 
cluster has been multiplied by a factor ~ 2. This is not to be linked with the dynamical 
evolution of the cluster, but with a much fainter background, i.e. a lower density. This 
last point demonstrates that extra caution should be taken when analyzing both the 
observations and the simulations of the morphology of the star clusters. 




Figure III. 21: Surface density profiles of the King9 model along the orbit B (at 
t = 20 Myr, red and t = 40 Myr, blue) compared to the isolated case (green). The 
vertical ticks mark the half-mass and half-light radii. The horizontal dotted lines denote 
the level of the local surface density of the surrounding galaxy. Obviously, such a level 
does not exist in the isolated case. 
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King?, t = 2Q Myr King?, L = 40 Myr 




Figure III. 22: Surface brightness of the three King models along the orbit B at the end of 
the compressive event (t = 20 Myr, left) and in the extensive mode (t = 40 Myr, right). 
The solid line circles mark the half-mass radii (considering all the stars) and the dotted 
circles show the half-light radii (from stars which brightness exceeds twice the level of 
the background). Each pixel covers 0.5 x 0.5 pc 2 . The colors range from white (surface 
density = 0) to black (central surface density of the cluster). 
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3.3 Conclusions 

As expected since the previous Sections, the fast transition from a compressive tidal 
mode to an extensive one leads to severe disruption and mass-loss of clusters yielding 
mildly peaked density profiles. With these first results, one can get a broader picture of 
the combined effects at stake at cluster scale: 

• As we have seen, a fast evolving tidal field impacts on the external layers of a 
cluster. This changes significantly the morphology, mass, luminosity profile and 
internal dynamics. 

• An initial mass function for the stars within the cluster also speeds up the evolution, 
via mass segregation (White 1977; Fleck et al. 2006). The more massive stars 
migrate to the center, and then weaken it through violent feedback (winds and 
SN blasts). (This point will be explored in the next months, with Nbody6++ and 
Titens.) 

These two effects combined act in the central regions and the external layers of a cluster 
by modifying its evolution. Therefore, a relation with the mass function of the clusters 
in the Antennae, and in merging galaxies in general, could be established. The question 
of "infant mortality" should be better understood too. 

Note also the importance of the binary/multiple fraction in the stellar population 
that is of prime importance when deriving velocity dispersions (see the recent works of 
Kouwenhoven Sz de Grijs 2008 and Gieles, Sana, & Portegies Zwart 2010), and thus the 
dynamical stage of a cluster (in particular, supervirialized or not). This aspect will also 
be considered in the next steps of our study. 
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Overview 

In this Chapter, we extend our conclusions about the 
Antennae galaxies to a broad range of configurations. 
Over the entire Chapter, the Antennae model presented 
before is considered as the reference for our parameter 
survey. We detail here the results initially presented in 
the Section 5 of Renaud et al. (2009). 



1 Need for a parameter survey 

During the first part of this thesis, we have focused on the Antennae galaxies. All the 
methods and diagnostics have been implemented with the view to compare the results 
obtained with direct observations of a real object. However, it seems natural to wonder 
if this pair of galaxies is representative of a larger set of configurations. 

Our results about the Antennae point out a strong correlation between the compres- 
sive tidal mode and the early life of star clusters. The next step of our study is thus 
to determine whether the same conclusions are applicable to a broad sample of inter- 
acting galaxies, and more precisely, which parameters are relevant for the growth of the 
compressive tides, and which are not. In the following Sections, we explore the role 
of several of these parameters by covering a significant volume of the entire parameter 
space, centered on those of the Antennae model (see Figure IV. 1). The main properties 
of the compressive tidal modes of these models are summarized in Table IV. 1. 

Our goal here is to compare the different models with the Antennae, that we consider 
and will refer to as the reference of this survey. 
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Table IV. 1: Parameter survey 



Model 


Key parameter Peak* 1 [%] Peak* 2 [%) 


7i* [Myr] 


Antennae (see Chapter II) 11.5 


11.5 


6.7 




Spin (see Section IV. 2) 






RP 


retrograde-prograde 7.5 


14.0 


7.3 


PR 


prograde-retrograde 7.5 


14.0 


7.3 


RR 


retrograde-retrograde 7.0 


18.5 


6.7 




Orbits (see Section IV. 3) 






OP 


high eccentricity 10.0 


9.5 


7.4 


OE 


low eccentricity 14.5 


8.5 


6.5 


OD 


distant encounter 3.5 


9.5 


7.3 


OC 


close encounter 15.5 




5.5 




Impact angle (see Section IV.4) 






IF 


face-on collision 39.5 


12.5 




IE 


edge-on collision 21.0 


8.5 


7.0 




Mass ratio (see Section IV. 5) 






M2 


2:1 mass ratio 8.5 


5.5 


6.6 


M3 


3:1 mass ratio 5.5 


6.5 


6.3 


M10 


10:1 mass ratio 1.5 


3.5 


5.5 




Progenitors model (see Section IV.6) 






KD 


Kuijken & Dubinski (1995) 4.5 


1.5 


3.3 


MD 


McMillan & Dehnen (2007) 10.0 


9.0 


2.6 



Note: Based on Tabic 3 of Renaud et al. (2009). 

* Mass fraction in compressive mode added by the first pericenter passage (compared to the 
isolation stage) . 

t Mass fraction in compressive mode added by the second major peak (merger phase). 
+ Exponential timescale of the lus distribution for the short time interval ATi. Note that 
this is not the same definition than in Renaud et al. (2009; recall Equation 11.25, page 82 
of the present document. 
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2 Spin 



As presented in Section 2.2 (page 8), the spin-orbit coupling of interacting galaxies is 
of paramount importance for their dynamical evolution, in term of timescale and creation 
of tails. For a galaxy B orbiting a galaxy A, the transition from prograde to retrograde 
goes through all the inclination angles % of the disk of B with respect to the orbital plane 
around A: for i — 0°, one has a prograde encounter, while i = 180° yields a retrograde 
system. In the following, we extend this definition to intermediate angles (i G [0°, 90°] and 
i G [90°, 180°]) so that a prograde encounter corresponds to an important perturbation 
of the disk by the tidal interaction, while a retrograde configuration only drives mild, 
transient perturbations. 
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Figure IV. 1: Schematic map of the original parameter survey presented in Renaud et al. 
(2009). The models in green have been added in the present document. Along each axis, 
a circles represents a model with a given value for the parameter considered. 



145 



Chap. IV 



Many interacting galaxies 



Toomre & Toomre (1972) presented a study on the influence of the spin-orbit coupling 
and concluded that it is a key factor in the creation of tidal tails. In this Section, we 
propose a similar work, but for self-consistent disks embedded in DM halos. Three 

PP | configurations are considered: in the PP (prograde-prograde) model, both galaxies have 
prograde encounters. This is the Antennae model presented in Chapter II. By flipping 

PR one of the disks, one gets the PR (prograde- retrograde) or the RP (retrograde-prograde) 
model. Note that the choice of flipping one disk or the other is not important, because of 
the high level of symmetry in the simulation. Finally, when the spins of both disks have 

RR | been inverted, one obtains the RR (retrograde-retrograde) model. The morphologies of 
the four models are shown in Figure IV. 2. 



2.1 Mass-loss 

The long tails existing in the Antennae galaxies are a strong hint that the encounter 
which created them was prograde for both disks. In this PP model, each disk ejects 
~ 30% of its mass in the form of tails, as already noted by Barnes (1988). Note that this 
value depends on the geometry of the merger and varies with, e.g. the pericenter distance. 
However, if one only changes the disk inclination (by steps of 180°), the mass-loss remains 
at the level of ~ 30% (prograde) or ~ (retrograde). 

A significant fraction of these ejecta will progressively fall back onto the nuclei, which 
helps the dissolution of the tails, long after the merger stage (~ 1 Gyr). This means that 
the tidal tails are not completely unbound to the central part of the progenitor galaxies. 
However, during the long range interaction period (i.e. between the first passage and the 
merger), the dynamical mass of the galaxies are strongly affected by the creation of the 
arms, which directly influences their trajectories. 

Figure IV. 3 shows the evolution of the distance between the two galaxies, for each 
model. First of all, no major difference is visible before the first pericenter passage 
(t < 0). Afterwards, the PR and RP models yield very comparable behaviors. For this 
reason, we will only consider the PR model in the following. Because of different amounts 
of mass-loss, the time of the second pericenter is shifted by ~ 50 Myr per tail created, i.e. 
it occurs at ~ 280 Myr for PP (two tails), ~ 330 Myr for PR (one tail) and ~ 380 Myr 
for RR (no tail). Keeping these dynamical data in mind, it is possible to make a link 
with the tidal field one expects to find in such systems. 



2.2 Compressive mode 

In Chapter II, we have seen that the major events in the evolution of the compressive 
tidal mode match the pericenter passages. Therefore, with the dynamics presented in 
the previous Section, one expects the mass fraction in compressive mode in the PR and 
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Figure IV. 2: Morphology of the PP, RP, PR and RR models at approximately the same 
dynamical step (maximum separation of the progenitors). The link between the prograde 
configuration and the creation of tidal tails is clearly visible. Here again, the red dots 
represent particles in compressive mode. 
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RR models to be peaked at the time of their respective passages. Indeed, as shown in 
Figure IV. 4, the previously mentioned shift of 50 Myr is visible between the three models. 

The main differences in this Figure correspond to the second pericenter passage, when, 
in addition to being shifted in time, the peaks have different amplitudes. Again this can 
be explained by a dynamical analysis. First, consider that 30% of the stellar mass of PP 
have been ejected to the tails during the first passage and thus, will not be affected when 
the nuclei interact again. Let Mo be the total initial stellar mass, and Mi, the mass that 
remains in the nuclei. In the PP case, one has Mi = (1 — 0.3)M = 0.7M . By reading 
the amplitude of the second peak, we see that ~ 15% of M , or (15%/0.7)Mi w 21% of 
Mi is in compressive mode. This fraction of the "bound" mass should be the same for 
the three models. For PR, only one tail is created, so Mi = (1 — 0.3/2)Mo = 0.85Mo. 
One expects 21% of Mi = 21% x 0.85M 17% of M to be in compressive mode at the 
second passage. This corresponds well to the value in Figure IV. 4. The same way, no 
tail is created in RR, so Mi M and the expected fraction is 21% of M , which also 
matches the value measured. 




Figure IV. 3: Evolution of the distance between the two galaxies, for the PP (black, as in 
Figure 11.15, page 76), the RP (green), the PR (red) and the RR (blue) models. RP and 
PR are very similar, due to the high level of symmetry of the progenitors. The differences 
in dynamical masses clearly affect the orbital period and separation. The two pericenter 
passages are marked with • and + respectively. 
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In addition, differences between the three models also exist between the two passages, 
corresponding to a re-arrangement of the matter within each galaxy according to proper 
internal dynamics. Elsewhere all models yield comparable behaviors, showing a quiet 
stage (intrinsic 3%), an oscillatory phase for the merger and finally a stable fraction with 
variations similar to the initial noise. All together, these points show that the major 
events in the compressive history of the mergers (peaks) are dominated by dynamical 
properties (mass and orbital period). 

Remains the analysis of the timescale of the compressive events. Do the spins of the 
galaxies influence the time a cluster-size element spends in compressive mode? We have 
seen that the lus and tt distributions depend on the time of integration, i.e. the time given 
to the particles to enter or exit the compressive regions (recall Section II. 3. 3, page 82). 
In the case of the Antennae model, this time has been arbitrarily defined and an order 
of magnitude for the timescale (~ 10 7 yr) has been obtained. For the comparison with 
the PR and RR configurations however, one requires a greater precision to highlight the 
influence of the spin on the timescale (which is expected to change slightly but remain 




Figure IV.4: Evolution of the mass fraction in compressive mode for the PP, PR and RR 
models. As before, t = corresponds to the first pericenter passage. A clear shift in time 
is visible between the models and can be simply explained by dynamical considerations. 
(Adaptation of Figure 10 of Renaud et al. 2009.) 
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at the same order of magnitude). To do so, the time of integration should be strictly 
the same for all three models. But as noted above, the dynamics of the systems strongly 
differs and thus, the times granted to the particles of the three models to enter and exit 
compressive zones are different. 

To balance this, we stop the integration at the time of the third peak in the merger 
phase (see the arrows on Figure IV.4). This gives a good equivalent stage in the dy- 
namical evolution of the merger. Note that it would have been preferable to stop later 
(e.g. at 1 Gyr, as in Chapter II) but it is not possible to precisely extract a dynamical 
stage from either the morphology (because the progenitors have merged) or the mass 
fraction in compressive mode. In this context, a comparison of the tt distributions is not 
relevant (because of its sensitivity to the time of integration), and one should focus on 
lus only. Figure IV. 5 shows that the differences between the three models represent less 
than 1%, and that all the properties of the distributions (double-exponential law, knee 
at ~ 50 Myr) are conserved when changing the spin. We note that the characteristic 
times Ti and T2 remain the same for all three models (see Table IV. 1). From Figure 11.29 




200 



Compressive mode interval At [Myr] 



Figure IV.5: lus distributions of the PP, PR and RR models, integrated from 
t = —100 Myr (isolated stage) to their respective merger phase. No clear trend is visible, 
which demonstrates that the orientation of the spin of the progenitors does not influence 
the time spent in compressive mode. 



150 



3 - Orbits 



Table IV. 2: Orbital survey 



Model 


e 


d [kpc] 


remarks 


Antennae 


0.96 


5.8 


reference model 


OP* 


1.10 


5.8 


almost parabolic 


OE 


0.80 


5.8 


low eccentricity 


OD 


0.96 


13.2 


distant encounter 


OC* 


1.50 


3.5 


close encounter 



Note: Table 4 of Renaud et al. (2009). 
* Orbital decay makes this orbit bound. 

' High eccentricity is required to avoid a central collision that 
would hide the effect of the initially non-zero pericenter distance. 



(page 93), we know that the tidal tails yield similar timescales for the lus distribution 
than the nuclei do (over ATi). Therefore, the presence of zero, one, or two tail(s) does 
not have a major influence on the lus distribution, as visible when changing the spin. 



3 Orbits 

In observed pairs of interacting galaxies, both the orbital period and the pericenter 
distance take heterogeneous values, because of (i) a higher density of galaxies at high 
redshift and (ii) few constrains on the geometry of the encounters. The exploration of 
the role of different orbital configurations of interacting pairs is thus important in our 
study of the tidal field. For the next part of our survey, we focus on the orbit of the 
progenitors, keeping all the others parameters (e.g. mass, spin, initial shape) as in the 
Antennae model. Therefore, the only changes made here involve the eccentricity e and 
the pericenter distance d. 

Because of dynamical friction, the centers of mass of the galaxies do not have a 
Keplerian motion and thus, it is difficult to describe their orbits accurately. For the sake 
of simplicity, the orbits investigated here are defined thanks to a Keplerian equivalent: 
from the isolated stage, we measure the eccentricity and the pericenter distance that the 
galaxies would have if replaced by point-masses, i.e. without orbital decay. Note that 
these quantities do not correctly fit the real values after the first passage (see Section A. 2, 
page 198) and are only used as an informational data to ease the comparisons between 
our models. Table IV.2 summarizes the properties of the five models presented here. 
(In total, 16 orbital configurations have been considered in this survey, but only five are 
presented here, for the sake of clarity.) 
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3.1 Eccentricity and distance 

Figure IV. 6 shows the morphology of the systems and the locations of their com- 
pressive modes when the galaxies have reached their maximum separation. In parallel, 
Figure IV. 7 displays the evolution of the distance between the two galaxies, and the 

OP | evolution of the mass fraction in compressive mode is shown in Figure IV. 8. The OP 
model develops very long tails as well as an important bridge. The compressive regions 
that exist in these tidal structures short after the first pericenter passage have already 

OE I disappeared in the snapshot shown in Figure IV. 6 (t ~ 450 Myr). The OE run shows a 




-50 50 -50 50 

x [kpc] 



Figure IV.6: Morphology of galaxies from our orbital survey, at approximately the same 
dynamical step (maximum separation of the progenitors). Note that the galaxies are 
not bound in the OC model. By contrast, OE yields a very short separation (see also 
Figure IV. 7). 
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very different configuration. Because of the low eccentricity of the orbit, a large fraction 
of the orbital angular momentum of the progenitors is transferee! to individual particles 
during the first passage (which dynamically heats the galaxies). As a consequence, the 
merger occurs just after the first encounter, forbidding a large separation of the galaxies. 
Thus, the transient tidal features (mainly in the tails and bridges) still exist when the 
progenitors start to merge (t ~ 150 Myr). 

As expected, the peaks of compressive modes match the pericenter passages, and 
the merger phase is also clearly marked. As for the Antennae, the mass fraction in 
compressive mode returns to a quiescent level (~ 3%) between the interaction events 
and remains stable for an arbitrarily long time after the merger. The Antennae and 
the OP models share a similar amount of mass in compressive mode at the first passage 
(t ~ 0) while OE yields a significantly higher value. By contrast with the spin survey 
(see the previous Section), the change already occurs at the first pericenter, i.e. at a 
stage when the disks are merely unchanged from their initial state. This demonstrates 
that the variation is not due to tidal structures, because they have not been created yet. 
In fact, in addition to the eccentricity, the relative velocity of the progenitors has a much 




200 400 600 800 1000 1200 

time [Myr] 

Figure IV. 7: Evolution of the distance between the two galaxies for the models of our 
orbital survey. Note that the vertical-axis is in logarithmic scale. The merger phase (for 
all models but OC) is visible as damped oscillations. 
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lower value for OE than in the Antennae and OP cases, suggesting that large compressive 
regions have time to develop all over the disks during the first passage. To conclude, the 
eccentricity is of paramount importance to set the characteristics of the orbits (period, 
friction) but does not have a strong direct impact on the tidal field itself. 

OD Figure IV.8 (bottom panel) clearly exhibits that the distant encounter (OD) provokes 
OC a smaller amount of matter to enter in compressive mode than the close encounter (OC). 
Furthermore, the second peak for OD (t ~ 700 Myr) corresponds to a second passage, 
much closer (~ 5 kpc) than the first one (~ 13 kpc at t — 0), and induces a higher mass 
fraction in compressive mode. This suggests that a close passage yields larger compressive 
regions than a more distant one. 



3.2 Flybys and mergers 

To investigate the last point further, we have conducted a secondary orbital survey, 
varying the pericenter distance only. Figure IV. 9 presents several measurements of the 



20 P~" r 




200 400 600 800 1000 1200 

t [Myr] 



Figure IV.8: Top: Effect of the eccentricity on the mass fraction in compressive mode. 
Bottom: Same but for the pericenter distance. (Adaptation of Figure 12 of Renaud et 
al. 2009.) 
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amplitude of the first peak as a function of the distance, covering very close passages 
and more remote encounters. In the later case, the progenitors undergo a mild inter- 
action and easily escape from their mutual gravitational attraction. Note that such a 
flyby configuration can be obtained with a high velocity, and yet close encounter. How- 
ever, this survey does not include high-velocity encounters, which limits its scope to 
slow (~ 200 km s _1 ) field-like galaxies (as opposed to clusters of galaxies having typical 
velocity dispersions of ~ 1000 km s _1 ). 

In the case of non-merging galaxies (right-hand side of Figure IV.9), the mass fraction 
in compressive mode remains at the same level than in isolation. That is, the interaction 
is too weak to induce compressive modes (other than intrinsic). When the distance 
decreases however, the fraction increases slowly before peaking sharply at d/Rd ~ 1, i.e. 
for distances comparable to the scalelength of the progenitors' disks. 
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Figure IV.9: Mass fraction in compressive mode at the first pericenter passage versus the 
impact parameter d, normalized to the disk scalelength R d (blue dots). For these points, 
the vertical line marks the limit between flyby and merger galaxies. The red symbols 
are measurements of maximum SFRs by di Matteo et al. (2007, their Figure 21). Filled 
triangles represent the mergers while the open ones denote the flybys. Both can exist for 
all d/Rd because of different relative velocities. (Adaptation of Figure 13 of Renaud et 
al. 2009, with data points from di Matteo et al. 2007.) 
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When linking the compressive modes to star formation, as in Chapter III, it is inter- 
esting to compare this trend with the results from simulations containing gas. di Matteo 
et al. (2007) performed a large survey of interacting galaxies with the aim to study the 
SFR. The red triangles on Figure IV. 9 represent the maximum value of the SFR they 
measured, for direct (prograde-prograde) encounters of spiral galaxies (that they called 
"L-L dir") only. The distance is normalized to their disk scalelength (R^ is taken equal 
to the parameter a* of their Miyamoto-Nagai models). 

For flybys (open triangles), the SFR tends to slightly decrease when the distance 
increases. Although in the case of mergers, these authors noted that the SFR increases 
with the distance. Over the range d/Rd ~ 2 — 4 however, the SFR seems to be fairly 
independent of the distance, at least for spiral-spiral mergers (filled triangles). It is 
important to note that the quantity considered here is the maximum SFR, i.e. non 
necessarily the value at the first passage. Therefore, the value of the SFR from di Matteo 
et al. (2007) can correspond to either a burst in star formation, or a slow increase due 
to, say, gas infall. Unfortunately, this does not allow to link star forming episodes with 
tidal events, as we do. 

With surveys like the one of di Matteo, making comparisons between the time evolu- 
tion of the SFR and the one of the compressive mode would certainly bring new lights 
on the exact role of the tidal field, with respect to gas dynamics like stripping or infall. 
The first steps of such a study have begun with our model of the Antennae and SPH and 
AMR runs (see Section III.2.5, page 122). 



4 Impact angle 

From the previous study of the influence of the spin and the pericenter distance, one 
knows that the importance of the compressive mode increases with (i) a high mass in- 
volved in the collision and (ii) a short range interaction. To go further on these points, 
two new configurations are presented in this Section 21 (see Figure IV. 10). The first one 

IF | (IF) is the face-on collision of two spiral galaxies, with disks perpendicular to the orbital 
plane, as in Section 1.1.6 (page 34). (Note that the progenitors are not raw exponential 
disks anymore but the usual models used in the Antennae reference.) In the second one 

IE | (IE), the disks are aligned with the orbital plane, so that they first collide through their 
edges. These two models should not be seen as statistically relevant for real configura- 
tions, but rather as limits for the relative inclination of the disks. 

In both cases, the pericenter distance is null and the disks overlap almost entirely 
at this time. That is, ~ 100% of the stellar mass sits in the configuration of a cen- 
tral collision, which is expected to rise large compressive regions (recall Figure IV. 9). 
Figure IV. 11 shows the evolution of the mass fraction in compressive mode of the two 

a These new models are not part of the survey published in Renaud et al. (2009). 
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Figure IV. 10: Schematic view of the configuration of the face-on (IF) and edge-on (IE) 
collisions. Red arrows indicates velocity vectors. 



models, compared to the Antennae. Interestingly, a small, minor peak is visible at t ~ 0, 
i.e. when the two nuclei overlap. Hence, the time of the first pericenter passage does not 
exactly correspond to the major peak in the mass fraction. However, the rapid increase 
occurs ~ 10 Myr after, when the nuclei are already separated by ~ 4 kpc. This matches 
the creation of very large compressive regions in the tidal tails created shortly (but not 
immediately) after the passage. As before, these compressive structures disappear within 
a few 10 7 yr, and new ones are created during the second passage and the merger phase. 

The amplitude of the peaks for both IF and IE models confirms that the compressive 
tidal mode scales with the amount of matter that goes "deep" into the potential wells. 
The lus timescales are comparable to those of the Antennae. However, the IF model is 
badly fitted with the usual double exponential profile, because of a very large number of 
transient features (representing ~ 40% of the stellar mass and lasting for ~ 5 Myr) in 
the tidal tails which expend at high velocity. Therefore, the lus statistics are biased by 
this short-period contribution. 



5 Mass ratio 

The cold dark matter theory suggests that the galaxy-size structures in the Universe | CDM 
form through the accretion of smaller aggregates, i.e. via repeated merger events (see e.g. 
Peebles 1982; Springel et al. 2005; Stewart et al. 2008 and references therein). The fine 
details of the evolution of such objects rely on the mass ratio between the components 
of the merger. At a smaller scale and when including baryonic matter, this parameter 
also seems to be important in the study of interacting pairs of galaxies (Naab & Burkert 
2003; Bournaud, Jog, & Combes 2005; Naab, Jesseit, & Burkert 2006; Johansson, Naab, 
& Burkert 2009). These works tell apart the major mergers (mass ratio greater than 3:1) 
and the minor mergers. However, it is important to precise which mass one refers to. 
Stewart (2009) noted that observational works often use the baryonic mass (which can 
be measured via light) while the theoreticians consider the baryons and the DM. In this 
study, we adopt the theoretical point of view and refer to the total mass (stars and DM) 
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Figure IV. 11: Evolution of the mass fraction in compressive mode for the face-on collision 
(IF) and the edge-on one (IE). 



when using the mass ratio. 

For this section of our survey, we have considered a major merger (2:1), a minor merger 
(10:1) and an intermediate one (3:1) in addition to our Antennae reference (1:1), and 
kept all the other parameters (spin, orbit, shape) as in the reference. However, changing 
the mass of the galaxies modifies the gravitational energy, while the kinetic energy of 
their particles is left unchanged. This inevitably leads to instability because a galaxy 
of lower mass but same velocity dispersion will dissolve. To avoid this, the size and the 
velocity have been scaled in addition to the mass, so that the density and the virial ratio 
are kept constant. This choice maintains the intrinsic mass fraction in compressive mode 
at ~ 3% of the disk mass. 



Box IV. 1 

When changing the mass of an iV-body system, one has to ensure that the 
internal dynamics are not affected. To keep the two-body heating low, the 
mass of individual particles should remain the same. Therefore, creating a low 
mass galaxy from an heavier one implies to reduce the number of particles. 
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Furthermore, when modifying the 


mass from Mi to M 2 , one can preserve 


the density p and the virial ratio v: 
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where £ and a represent the length and velocity scales. This leads to 
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The morphology of the galaxies is shown in Figure IV. 12. On the one hand, for the 
major and intermediate mergers (M2 and M3), the morphology and the position of the M2 
compressive modes in the southern galaxy (which has not been changed) remain very M3 
similar to those of the Antennae. In the minor merger case (M10) however, the small M10 
intruder clearly has a too weak influence on the main disk which only exhibits a faint 
tidal tail. On the other hand, the modification of the aspect of the counterpart strongly 
depends on its relative mass: the low mass progenitors undergo mass loss via tidal 
stripping. In the same time, they induce radial perturbations in the main disk, which 
lead to the creation of transient features like spiral arms and even a bar. We note that 
such structures generally yield cores in the gravitational potential and thus compressive 
regions, as visible in the lower-right panel of Figure IV. 12. As a consequence, the mass 
fraction in compressive mode is expected to remain fairly low (close to isolation level) 
in the satellite galaxy but higher in the major progenitor, than in the nuclei of 1:1 
encounters. 

This idea is confirmed by Figure IV. 13. As always, the first passage is clearly marked 
by a peak. For the minor mergers (M3 and M10), the decrease of the mass fraction 
between the two passages is replaced by a slow increase, due to the creation of the 
transient compressive structures (spiral and bar) mentioned above. In the M3 model, 
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Figure IV. 12: Morphology of the Antennae, M2, M3 and MIO models at approximately 
the same dynamical step (maximum separation of the progenitors). The minor merger 
leaves the main disk (south) merely unaffected. However, transient spiral structures are 
created during the interaction. 
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they still exist when the second passage occurs (t ~ 425 Myr), while they have time to 
slowly vanish in M10 until the merger (t ~ 900 Myr). In this case, the maximum amount 
of matter in compressive mode is reached between two interaction events. However, it 
is still possible to tell it apart, as it does not correspond to a peaked, rapid change. In 
terms of timescale, this can be compared to the transient compressive modes found in 
the tidal tails and the bridges of the Antennae (recall Figure 11.29, page 93). 

In fine, minor mergers induce many short-lived compressive regions, leading to a 
sharper lus distribution. As a result, the exponential timescale T\ is slightly shorter for 
the minor mergers than for the major ones. By contrast with the orbital survey, these 
conclusions derive from internal dynamics and are not biased by the orbital period. 
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Figure IV. 13: Evolution of the mass fraction in compressive mode for the major (An- 
tennae, M2), intermediate (M3) and minor (M10) mergers of our survey. Despite very 
different dynamics, the relation between the interaction events and the peaks is still clear. 
(Adaptation of Figure 14 of Renaud et al. 2009.) 
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6 Further compound models 

Up to now, we have always used an exponential disk, an Hernquist (1990) bulge and 
an isothermal halo for our progenitors. However, the halo is suspected to influence the 
characteristics of the tidal field, because it dominates the gravitational potential of the 
galaxy. To verify this hypothesis, the last step of this parameter survey sets up different 
compound galaxy models. However, it is not possible to simply modify the parameters of 
the halo, leaving the other two components unchanged. In composite models, the global 
potential is built by the sum of the three parts. Therefore, changing one requires the 
adjustment of the others, in order to conserve the stability and coherence of the galaxy. 
That is why the disk and the bulge have to be modified, too. To do so, equilibrium 
distribution functions from the literature 5 have been used to create two new composite 
models, simply keeping a circular velocity of ~ 220 km s -1 at solar radius. 

The first one (noted KD) is a Kuijken & Dubinski (1995, model MW-A) model of a 
Milky Way-like galaxy. It gathers a 3D version of a Shu (1969) disk, a King (1962) bulge 
and a lowered Evans (1993; see also Kuijken & Dubinski 1994) DM halo. The iV-body 
realization of this model have been done with the mkkd95 tool of the Nemo package. KD 
has slightly more massive halo (xl.3) and disk (xl.5) than the reference model. The 
scaling is so that the truncation radius of the disk is the same than in the reference. 

These changes in the shape of the galaxy imply modifications of the potential and thus 
of the tidal profile too. Figure IV. 14 plots the maximum eigenvalue of the tidal tensor 
as a function of the radius, for the galaxies taken in isolation. The intrinsic central 
compressive region exists for all models and is clearly marked on this plot. We note that 
its diameter and tidal strength vary significantly from one model to the other. That is, 
one may expect different results from the (by-now) traditional study of the tidal field of 
the merger. These differences are not visible in the location of the compressive regions, 
as shown in Figure IV. 15: the major areas of compressive tides (nuclei, bridges, tails) 
are retrieved in a comparable way than in the reference model. 

Figure IV. 16 displays the evolution of the distance between the two progenitors, mark- 
ing the usual steps: first passage, separation, second passage and merger. The three 
configurations lead to a very similar evolution, especially between the reference and KD. 
In parallel, Figure IV. 17 shows the mass fraction in compressive mode. As expected, the 
intrinsic fraction (i.e. from isolation stage) is different for KD from the Antennae case. 
This is explained by (i) a smaller central compressive region and (ii) a different density 
profile (or in other words, a different mass enclosed inside the compressive limit). Only 
~ 0.6% of the disk's mass is intrinsically in compressive mode for KD, instead of ~ 3% 
for the Antennae model. However, we note a very good correlation between the two 

b The details on the construction of such distribution functions are out of the scope of the present 
document. The interested reader will find their derivations in the references presented below. 

c Note that the conversion factors from TV-body to physical units are slightly different than previously: 
3.6 kpc, 5.5 x 10 10 M Q and 14.2 Myr. 
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Figure IV. 14: Absolute value of the maximum eigenvalue of the tidal tensor of the three 
compound models (Antennae, KD and MD, in isolation), as a function of radius. Dashed 
curves represent the negative (i.e. compressive) part of the profiles. 
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Figure IV. 15: Morphology of the KD and MD models at approximately the same dy- 
namical step (maximum separation of the progenitors). 
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models in the evolution of this fraction through the merger. In KD, the peaks are better 
separated in time which suggests a shorter timescale for the compressive mode, and a 
more efficient orbital mixing of the particles. This is confirmed when investigating the lus 
distribution: a timescale of ~ 3.3 Myr is measured, instead of ~ 7 Myr (see Table IV. 1). 
Here again, the double exponential distribution is very well retrieved. 

MD The last model of our survey is taken from McMillan & Dehnen (2007, hereafter MD). 
NFW This time, the halo is a massive NFW profile (Navarro, Frenk, & White 1997). It is 
combined with an Hernquist (1990) bulge and an exponential disk. The iV-body model d 
is created thanks to the mkgalaxy tool within Nemo. The major difference between this 
model and the reference is the extension of its halo: it shows a characteristic radius of 
12 kpc and is 24 times more massive than the one of the Antennae model. 

As for KD, the intrinsic mass fraction in compressive mode (see Figure IV. 17) is quite 
low (~ 0.7%). In addition to the central compressive region, number of small, transient 
compressive zones exist all over the disk, because of a peculiar local distribution of the 



d As before, the conversion factor from A^-body to physical units are different than in the Antennae 
case: 2.8 kpc, 6.0 x 10 10 M Q and 8.8 Myr. 
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Figure IV.16: Evolution of the distance between the two galaxies for the KD and MD 
models. Except a faster evolution of the MD model, the global behaviors are comparable 
to those of the Antennae reference. 
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mass. Note that these features are smoothed out when averaging the tidal profile per 
radial bins, as in Figure IV. 14. The massive halo and the important dynamical friction 
reduce the orbital period by ~ 130 Myr (Figure IV. 16), which forbids the mass fraction 
in compressive mode to return to its isolation level between the two passages, like in the 
M3 model presented in the previous Section. 

We note that the NFW halo of a progenitor extends well beyond the limits of the 
disk, so that the entire stellar structure (initially the disk, and then the nucleus and the 
tails) sits in the inner part of the NFW profile (see Section B.1.5, page 208). There, the 
density profile of the halo only (p oc r _1 ) yields less extensive tides than in the outer 
regions (p oc r~ 3 ), meaning that the additional effect of the disk and the bulge would 
more easily create compressive zones than in the external part of the halo. 

As in the KD model, the timescale of compressive modes is shorter than the one of 
the Antennae (~ 2.6 Myr) but again, the double exponential law is well reproduced. 

These two models confirm that the importance of compressive modes in term of mass 
fraction depends on the distribution of the matter in the progenitors. However, its evo- 
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Figure IV. 17: Evolution of the mass fraction in compressive mode for the three compound 
models. Differences in the intrinsic fraction (t < 0) and in the time lapse between the 
peaks exist, but the overall behavior is well conserved. 
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lution seems fairly independent of the internal structure of the galaxies, as the behavior 
of the tidal field follows the chain of events in the interaction. To go further in the 
investigation of the role of the halo only, Chapter V presents an analytical derivation of 
the tidal fields of typical halos, derived from cosmological simulations. 



7 Conclusions of the survey 

The survey presented in this Chapter allows an exploration of the role of the major 
parameters of a galaxy merger. The overall conclusion has three folds: 

• the formation of compressive regions is systematically triggered by major inter- 
action events (passages, merger phase). Flybys do not yield a compressive tidal 
field. 

• the mass fraction in compressive mode is strongly anti-correlated with the peri- 
center distance: a close encounter induces a large amount of matter to sit in such 
modes while a distant one has a milder influence. 

• the duration of these modes depends on the dynamics of the merger and of the 
individual galaxies, but always yields the same order of magnitude (10 7 yr). 

Finally, the conclusions drawn in Chapter II about the special case of the Antennae 
galaxies can be extended to a broad range of parameters. The orders of magnitudes 
derived apply to many (if not all) configurations and, more detailed descriptions follow 
the "rules of thumb" mentioned above. 
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Overview ^ 

In this Chapter, we widen our field of view to larger 
scales. Using the results of cosmological studies, we 
investigate the influence of the density profile of galactic 
halos on the character of the tidal field. This work is 
summarized in Renaud, Boily, & Theis (2010). 

) 



1 Dark matter halos 

1.1 Missing mass 

The early measurements of velocity dispersions of stars (Oort 1932) and galaxies in 
clusters (Zwicky 1933) revealed that the corresponding dynamical mass is ~ 10 — 100 
greater than what is derived from the light: the "luminous" mass can not explain the dy- 
namics of astrophysical objects (in the framework of the Newtonian physics). Therefore, 
it has been naturally suggested that this "missing mass" is not luminous, or in other 
words that DM exists in the Universe. Many decades after, the nature of this matter 
remains unknown and puzzling for physicists (see Trimble 1987; Carr 1994 for reviews). 
It is supposed that an important part of its mass would be in a form of non-baryonic mat- 
ter 3 , with the rest being non- luminous baryons (black holes, brown dwarf stars, planets). 
Hereafter, we take the shortcut of calling "baryonic", all the non dark matter. 

Detailed studies on individual galaxies revealed that their luminous mass (stars and 
gas) is embedded in a DM halo (see e.g. Freeman 1970). However, the distribution of this 
material is not well-constrained by observations, de facto. That is the reason why one has 

a These candidate particles are called WIMPs (weakly interacting massive particles) because they 
should have a very small cross-section to avoid interaction with the "classical" baryons that would make 
them easier to detect. 
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to assume a shape for the density profile of the DM halo, when modeling galaxies (isolated 
or interacting). In the previous Chapters, we have considered isothermal (Chapter II), 
Evans and NFW (Section IV. 6, page 162) profiles, via distribution functions found in the 
literature. Other profiles can be studied analytically and/or numerically in principle, but 
the distribution functions of composite models (i.e. containing a stellar disk and bulge 
in addition to the DM halo) are generally not available. Furthermore, the equilibrium 
conditions do not couple a distribution of the baryonic matter for every halo and thus, 
a disk embedded in a given halo cannot (in principle) be exported to another one. That 
is, measurements of the effect of the halo on a property of a galactic system (like the 
tidal field, as discussed in Section IV. 6, page 162) is often hidden by (if not buried in) 
the effect of different disks and bulges. 

In this Chapter, we circumvent this issue by investigating the properties of the DM 
halo only, with an analytical method. The goal of this approach is to highlight the 
differences existing between the profiles, in particular for the character of the tidal field 
(extensive or compressive). 

1.2 Cosmological simulations, cores and cusps 

The bottom-up scenario suggests that structures (and among them, the galaxies) form 
through the repeated accretion of clumps of matter. Following this idea, cosmological 
simulations have been able to reproduce the organization of galaxies in groups and clus- 
ters, as observed (see e.g. Tegmark et al. 2004; Springel et al. 2005). With such 
numerical studies, it is possible to retrieve the density profiles of the halos, and possibly 
derive an analytical definition. That is how Navarro, Frenk, & White (1997, hereafter 
NFW) fitted their simulated DM halos (once they have reached equilibrium) with the 
formula 



where pn,o is a characteristic density and rN,o a scale radius. This profile provides a good 
approximation on a wide range of radii and masses (as noted by e.g. Navarro, Frenk, & 
White 1997; Navarro et al. 2004), a strong hint that it is universal. Comparable results 
M99 I have been obtained by Moore et al. (1999, hereafter M99) who found 



Pn,o 



(V.l) 




2 ' 



Pm{t) 



Pm,o 



(V.2) 




Let's first focus on the logarithmic slope of these profiles, defined as 



dln(p) 
dln(r) 



(V.3) 
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Figure V.l plots this slope as a function of r. Both profiles yield (3 = 3 in their outer 
regions (i.e. for r ^> r ), and a finite value (1.0 and 1.5, respectively) for r <C r . 
However, more recent high resolution runs noted a small, yet systematic deviation of 
simulations from the analytical description (Navarro et al. 2004; Boily & Athanassoula 
2006; Hansen & Moore 2006). Indeed, there is no sign of convergence toward a finite 
value of [3 in the inner region, as already suggested by the data of Navarro et al. (2004, 
dashes on Figure V.l). In other words, several models define a density profile with a 
central cusp denoted by (3{r = 0) > 1, while simulations tend to predict {3(r = 0) = 0, 
i.e. cored halos (see de Blok 2010). 



Box V.l 

What distinguishes a core from a cusp? In spherical symmetry, a profile 
yields a central core if the slope of its potential at r = is null. In mathe- 
matical terms, a cored distribution verifies 

Km 1=0 

r-s>o V dr 



lim ( r 2 / x 2 p(x) dx 1 = 
Jo / 



r 

x 2 p(x) dx r ~° r k with k > 2. (V.4) 

Thanks to the Maclaurin expansion (i.e. Taylor series at r = 0) of the inte- 
grand, one has 

p(r) ~ r k 6 
/3(r = 0) < 1. (V.5) 

Therefore, all density profiles with a logarithmic slope /3 strictly smaller than 
unity correspond to cored potential. All the other cases yield a cusp. 



A more sophisticated profile 5 , 
halos (Navarro et al. 2004): 



called Einasto (1965), better matches the simulated 



Pe{t) = PE,oexp <^ - 



a 



rE,o 



(V.6) 



This profiles yields a logarithmic slope which is a power-law of radius c : (3(r) oc r a , 
so that the central part exhibits a core. When applied to the simulations, the best fit 

b The reference is difficult to find but we refer the reader to Cardone, Piedipalumbo, & Tortora (2005) 
who proposed a detailed investigation of the properties of the Einasto density profile. 

c Such a function of the volume density is equivalent to the Sersic (1968) law on surface density (see 
Merritt et al. 2005 for a discussion on this similarity). 
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is found for the index a = 0.17 (Navarro et al. 2010), a value that we adopt hereafter. 
Note that Stadel et al. (2009) found a = 0.15. The conclusions presented below are 
not affected by such a change and thus, we leave the discussion on the value of a to the 
cosmologists. 

Unfortunately, the difference between cored and cuspy profiles raises at small radii only 
and it is still difficult (with current resolutions) to distinguish them. In Section V.2, we 
show that their tidal fields however, are very different. 




Figure V.l: Logarithmic slope of the density profiles of galactic halos (r = 20 h~ l kpc) 
from cosmological simulations (Navarro et al. 2004, dashes), compared to the NFW 
(blue), M99 (purple) and Einasto (red) analytical curves. The slope of the Hernquist 
(green) profile is also plotted, for reference only. (Simulation data points from Navarro 
et al. 2004.) 
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1.3 7-family 

Before these simulations, Dehnen (1993) and Tremaine et al. (1994) analytically 
defined a family of density profiles as 

(r) = QzjLH M (v.7) 

47T r 7 (r + r 7i0 ) 4 7 

where r 7)0 is a scale radius and 7 G [0, 3[ a parameter. Note that for 7 = 1 and 7 = 2, 
one retrieves the Hernquist (1990) and Jaffe (1983) profiles, respectively. This definition 
has been first used to describe elliptical galaxies and bulges, but has been extended to 
DM halos later on (see e.g. Merritt & Fridman 1996; Jeon, Kim, & Ann 2009). Its main 
advantage for our topic is that it yields a wide range of logarithmic slopes: {3(r = 0) = 7, 
from cored configurations (7 < 1) to strong cusps. 



1.4 Making comparisons 

In order to compare the properties of several DM halos profiles, one has to ensure 
that they are properly scaled. Indeed, all the analytical forms use scale-free parameters 
(po, To, M) that do not always represent the same physical quantity. It is therefore of 
prime importance to establish a relation between them, before making any comparisons. 
Being halos of galaxies, all the profiles should have a comparable impact on the dynamics 
of the embedded disk. For this reason, we define a reference radius r s containing 1/5 
of the total mass (see e.g. the model of the Milky Way in Binney & Tremaine 2008), 
so that the circular rotation speed of the baryonic material at this point remains fairly 
unchanged. To do so, the integrated mass is computed for each profile, which gives 
the total mass when evaluated at an infinite radius. The mathematical derivations are 
detailed in Section B.l (page 205). Then the equation of the mass profile is solved to 
find the relation between tq and r s . 

In the case of Einasto, one finds (numerically) that r E ~ 0.36 r s . Then, p E can be 
determined: 



2 (3/a) a (l-3/a) / 2 \ 

" M 4, (0.36 r s rT ( r XP (") <V - 8) 
M 

0.14 — . (V.9) 



' s 



Repeating this exercise for the 7-family, one gets 

, fl = r s [b 1 ^ - 1] . (V.10) 



Thanks to Equation B.23 (page 208) and Equation B.37 (page 210), Figure V.2 shows 
the mass profiles for Einasto (red) and Hernquist (7 = 1, green). 
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Figure V.2: Integrated mass of the NFW, Einasto and Hernquist profiles as a function 
of the radius, normalized to r s so that M(r B ) = M/5. 



The case of NFW is more problematic, as its mass (see Equation B.30, page 209) 
integrates to infinity, which means that a total mass cannot be defined. To overcome 
this, one must recall that the NFW and Einasto profiles are comparable over a certain 
range of radii. Therefore, it is possible to set the scale of NFW via the one of Einasto. 
Figure V.2 plots the best fit of the Einasto mass profile with the analytical form of a 
NFW, which gives 



and 



r Ni0 « 0.34 r s 
~ 0.9 r Ej0 , 

p N ,o ~ 4.5 p E ,o- 



(V.ll) 



(V.12) 



The divergence for r>l leading to an infinite mass is clearly visible but one also notices 
that the two profiles are remarkably comparable for smaller radii. 

Figure V.3 displays the relative difference between the three profiles: for Einasto and 
NFW, it reaches the level of 30% for r w 40 r s only, i.e. for rj 0.61 M. Therefore, it is 
still difficult to distinguish between both models from simulations, which are limited in 
resolution. However, these small differences strongly impact on the associated tidal field. 
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Figure V.3: Relative difference between the integrated masses of the NFW and the 
Einasto profiles (blue), and between Hernquist and Einasto (green). 



2 Tidal field 



2.1 In isolation 



With scaling relations between the models, it is possible to go further in their com- 
parison and exhibit their tidal field. The analytical derivations are given in Section B.l 
(page 205). For simplicity, we will use M = 1 and r s = 1 in the following. 

First, it is noticeable that compressive tides only exist in the core of a potential, 
therefore, in the region where the logarithmic slope of the density profile is smaller 
than unity (recall Box V.l). With Figure V.l, it immediately appears that NFW, M99, 
Hernquist and the isothermal profiles do not have compressive tidal modes, while Einasto 
does in its inner part. 

As a confirmation, Figure V.4 displays the value T xx of the tidal tensor of these 
profiles, as a function of radius. As expected, only the Einasto profile yields compressive 



175 



Chap. V 



Tidal field and dark matter halos 



1 5 



T 



T 



T 




NFW 
Einasto 
Hernquist 



0.001 



0.010 



0.100 
r/r s 



1.000 



10.000 



Figure V.4: Tidal profile of the NFW, Einasto and Hernquist profiles. The horizontal 
dashed line marks the transition between extensive (> 0) and compressive (< 0) tides. 

tides at small radii. Note that even if this central part is very small d , the overall shape of 
the tidal profile strongly differs from the one of e.g. NFW, though their density profiles 
are comparable. 

The same way, Figure V.5 shows the tidal profile of several members of the 7-family. 
Compressive tides only exist if 7 = /3(r = 0) < 1. In these cases (corresponding to cored 
potentials), the compressive central region covers a radius of r 7j0 (1 — 7)/2. In all the 
other cases, the tides are extensive everywhere. Again, a broad range of configurations 
for the tidal field is accessible from a restricted amount of density profiles. 

2.2 In symmetrical pairs 

The analysis of spherically symmetric mass distributions reveals that the group of 
profiles commonly used for modeling the DM halos spawns very different tidal fields. 
As presented in the previous Chapters, the importance of the tides is highlighted when 

d In the typical case of the Milky Way, one has r s ~ 10 kpc. Therefore the compressive region of the 
equivalent Einasto profile would cover ~ 100 pc only. 
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Figure V.5: Same as Figure V.4, but for some models of the 7-family. 



galaxies interact. Therefore, this Section repeats the analysis made in the previous one, 
but when two profiles are involved. Unfortunately, the symmetry is lost and thus, a 
purely analytical derivation of the tidal field becomes (very) challenging 6 . The first 
possible solution is to perform iV-body simulations. However, such a method introduces 
noise in the rendition of the potential, which propagates in the computation of the tidal 
tensor. Furthermore, the Einasto and NFW profiles would be hardly distinguishable 
at typical resolution, and thus their differences would be buried in the iV-body noise. 
Finally, it is not always easy to construct distribution functions for all theoretical models, 
as mentioned before. 

The alternative adopted here consists in a semi-analytical method: the tidal tensor of 
each individual halo is calculated analytically (i.e. by hand) and then implemented as a 
function in a numerical code. Using the parameters of the global system (masses, sizes 
and positions of the halos), the code sums the individual tidal tensors at every point in 
space to get a global tensor and finally set it in diagonal form. Therefore the tidal field 
is computed for a static configuration (see Section V.3.2 for a discussion on this point 



c The derivation of the tidal tensor can still be done using the superposition principle for the potential, 
and the distributivity of its second derivative over the addition. However, the global tensor yields an 
azimuthal dependence that would be averaged if considering the radial profile only. 
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and the limits of this approach). We focus on configurations close to those of the first 
pericenter passage of interacting galaxies: the overlap of two identical potentials, one 
shifted by a distance d with respect to the other. 

Figure V.6 plots the maximum eigenvalue A max of the total tensor of the 7 profiles, 
with red denoting fully compressive tides and blue, extensive regions. In the immediate 
vicinity of the centers (white plus signs), the tidal field is as in isolation, i.e. a compressive 
mode only exists for 7 < 1. At larger distances however, non- intrinsic red zones appear, 
due to the overlap of the two halos. Interestingly, these zones are visible up to 7 ~ 2, while 
isolated models are totally extensive for 7 > 1. This indicates that the overlap of the 
potentials changes the very character of the field. Note that these regions are not situated 




Figure V.6: Map of the maximum eigenvalue of the tidal tensor for two identical 7 
profiles separated by a distance d = 0.4 r s and centered on the plus signs. Dotted circles 
mark the radii r s . The number in the top-right corner of each panel represents the value 
of 7. Black is used to enhance the contrast and denotes extremely compressive areas 
(A max < -!)• (Figure 1 of Renaud, Boily, & Theis 2010.) 
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Figure V.7: Same as Figure V.6 but for two NFW (left) and two Einasto (right) profiles. 
Note that the amplitude of the color bar has changed, to enhance the contrast. Clearly, 
the character of the tidal field (compressive or extensive) strongly depends on the details 
of the profile chosen. (Adaptation of Figure 2 of Renaud, Boily, & Theis 2010.) 



between the wells but on the perpendicular axis. Therefore, these compressive modes 
may impact on the (possibly) embedded baryonic matter in different ways, depending on 
its alignment with respect to the orbit of the merger. 

In Figure V.7 the experiment is repeated but with two NFW and two Einasto profiles. 
As before, the intrinsic compressive region of the Einasto model is found at the center of 
the wells. Note that its diameter (~ 10~ 2 r s , recall Figure V.4) corresponds to ~ 0.5 mm 
on the Figure. Therefore, it is hidden behind the plus sign f . In addition, the patterns 
found in Figure V.6 for e.g. 7 = 1.25 are approximatively reproduced, in sharp contrast 
with the NFW case that does not yield any compressive region. 

To conclude, the differences in the intrinsic tidal fields of Einasto and NFW already 
noted in Figure V.4 are strengthened when another, identical profile is involved: sec- 
ondary compressive regions exist between the Einasto wells while the entire space remains 
in extensive mode with NFW halos. 



f In the electronic version of the manuscript, the reader can zoom-in on the sign and discover a 
small, yet visible red and black pattern. With the highest zoom factor, it is even possible to detect an 
asymmetry: the extremely compressive part (black) is preferentially found opposite to the second well, 
while the red zone faces it. 
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2.3 Other configurations 

The results obtained in the previous Section can be extended to other distances d. 
First, we have checked that NFW profiles never give out compressive tidal modes, for any 
separation of the wells. Secondly, the upper row of Figure V.8 plots the map of A max for 
two Einasto halos closer (left panel) and more distant (right) than in the previous case. 
The size and amplitude of the non-central compressive regions show a clear dependence 
with d, as expected according to the results of Section IV. 3. 2 (page 154). 

To pursue the investigation of other configurations, the mass ratios /i = Mi/M 2 can 
also be changed, by setting a different weight for the two tidal tensors: T = Ti + \i T 2 . 
(Note that this method allows us to scale the NFW halo although it has an infinite 
mass.) With the view of modeling galactic halos, the radial scale (ro) is adjusted so that 
the average density is kept constant (recall Equation IV. 3, page 159). As expected, the 
NFW halos remain entirely in extensive mode. The bottom-row of Figure V.8 shows the 
results for Einasto profiles with the mass ratios 3:1 (left panel) and 10:1 (right panel). 
The symmetry of the compressive regions is obviously broken but the overall topology 
is well conserved. These regions stand closer to the lightest halo (the upper one on the 
Figure) than to the major one because the extensive modes are weaker there, and can 
easily be compensated or even overcome. 

The amplitude of both the compressive and the extensive modes is comparable in all 
three cases presented (1:1, 3:1 and 10:1). This confirms that the strength of the tidal 
field mainly depends on the distance between the structures and is fairly invariant with 
the mass ratio. 

The situation depicted above corresponds to well-defined events in cosmological simu- 
lations when for instance, a sub-halo falls onto a main object. In that case, the compres- 
sive modes that develop around its trajectory may influence the dynamics of the material 
already present in these regions. 



3 Validity and limits of the method 

3.1 Role of the baryons 

In this Chapter, we have focused on configurations applicable to interacting spiral 
galaxies for which the separation is a fraction of the diameter of the disks. At these 
scales, a sole description of DM is not enough to address the question pertaining to the 
dynamics of the system. Therefore, baryonic physics should be included, the same way 
it is done in the previous Chapters, by taking the bulges and the disks into account. 
However, the purpose of this study is to show that the halos alone already influence the 
character of the tidal field over ~ 10 kpc scales. 
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Figure V.8: Same as previous but for Einasto profiles separated by d = 0.2 r s (top-left) 
and d = 0.6 r s (top-right). In the bottom row, the distance remains d = 0.4(r Si i + r s 2 )/2 
but with a mass ratio of 3:1 (bottom-left) and 10:1 (bottom-right). In these cases, the 
lighter halo is the upper one. The global center of mass is always set at the origin. 
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The results presented demonstrate that, on the one hand the choice of the shape of the 
halo has a limited impact on the overall density profile but, on the other hand it strongly 
affects the properties of the tidal field. When putting these facts in the context of star 
formation (recall Chapter III), the NFW profile does not favor the formation of clusters 
or TDGs, though the comparable Einasto profile does. Note that we do not claim that 
a given profile would prevent star formation, but simply that its tidal field creates a less 
favorable environment with respect to another DM distribution. The additional baryonic 
physics would set the fine details of these processes, on the base of the conclusions drawn 
here. 



3.2 Static versus live 

This semi-analytic approach uses the superposition of two spherically symmetric halos 
separated by a distance of the same order of magnitude than their characteristic scale. 
Because this situation is comparable to what happens in interacting galaxies, one may 
ask whether the deformation of one halo due to the tidal field of the other can really 
be neglected. Indeed, the phenomenon responsible for the creation of the tidal tails also 
influences the morphology of the DM g . 

One of the assumptions of this study is that the static configurations presented resem- 
ble the one of a first pericenter passage, thus before the halos are strongly affected by the 
interaction. To verify this last point, Figure V.9 compares the maps of A max of a "real" 
live merger to the one from a static equivalent. In a first experiment, two Hernquist 
potentials are rendered with 5 x 10 5 particles each, and placed on an elliptical orbit. The 
equations of motions are solved by gyrf alcON (see Chapter I) up to the first pericenter 
passage, when the tidal field is computed (panel a). For the second experiment, the two 
halos are replaced by their initial mass profiles (iV-body again) at their positions at the 
first passage, enforcing de facto the spherical symmetry (panel b). The two maps are 
hardly distinguishable, which confirms that the halos are not morphologically evolved 
before the first passage. On panel c, the result of our semi-analytical approach for the 
same configuration is shown, which helps to appreciate the level of numerical noise in 
the previous two cases. 

3.3 Triaxial halos 

The DM halos extracted from cosmological simulations are not perfectly spherical 
objects but triaxial, ranging from oblate to prolate (see e.g. Davis et al. 1985; Hayashi, 
Navarro, & Springel 2007 and references therein). This asymmetry becomes important 
when considering a disk in rotation inside the halo: its plane should lie close to the one 

B in a milder way though, because the velocity dispersion of the halo is more important than the one 
of the disk. 
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Figure V.9: Map of A max for two live (a), static (b) and analytical (c) Hernquist potentials. 
The three experiments yields comparable results (except the numerical noise), which 
demonstrates that the halos are not morphologically affected by the interaction prior to 
their first passage. Note that the difference visible in the central part of the wells comes 
from a mathematical discontinuity at r = (see Equation B.20, page 207) that cannot 
be reproduced numerically. 
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Figure V.10: Map of A max when an elongated structure is anti-aligned (left) and aligned 
(right) with another potential. The patterns of compressive zones remain comparable to 
those given out by spherical models. 



defined by the axes of the halo. Not taking this into account would false the interpretation 
of the velocity curves of such disks and, in the end may even bias the conclusions on the 
inner slope of the DM profiles (see Hayashi & Navarro 2006 on the cusp-core problem). 

For a preliminary study like the one presented in this Chapter, we have restricted the 
investigation to spherically symmetric distributions, i.e. described by a monopole term. 
It is possible however to extrapolate to quadrupoles, or "bared" systems by combining 
two identical potentials A and B, separated by a distance e <C r s , thus forming an 
elongated structure. This new object has been used as one of the halos involved in our 
previous configurations, the other being a third, spherical potential C. To keep the mass 
ratio between the halos equal to unity, the A and B profiles have been downscaled by a 
factor two in mass (and 2 1 / 3 in size, recall Box IV. 1, page 158). The distance between 
the barycenters of AB and C is set to d — 0.4 r s ^> e, as before. 

Figure V.10 displays the usual map when the bar is anti-aligned and aligned with C. 
For both orientations, the resulting map is equivalent to the one of two pairs of potentials 
with respective separations e and d. The morphology of the compressive regions slightly 
changes with the orientation of the bar. However, the large scale structures remain fairly 
unaffected by the triaxiality. Separations up to e ~ d have also been explored to confirm 
the robustness of these points. In conclusion, the presence of a triaxial halo modifies the 
size and shape of the tidal regions (extensive and compressive) but does not influence 
the very character of the field. 
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4 Modified dynamics 

At the beginning of this Chapter, we have introduced to concept of DM to explain 
velocities higher than expected from the luminous mass in the context of Newtonian 
dynamics. However, instead of changing the mass, it is also possible to change the laws of 
physics. This is the line pursued by the modified Newtonian dynamics theory (Milgrom | MOND 
1983; Bekenstein & Milgrom 1984). It proposes that the intensity of the gravitational 
force F decreases as the inverse of the distance (V -1 ) in the weak field regime, while it 
stays Newtonian (r -2 ) elsewhere: 

F = m a ( —) a with li(x) = < ^ ^ \ (V.13) 
\a J { x if x < 1, ' 

where a is the acceleration and ao ~ 1.2 x 10~ 10 m s -2 the critical acceleration (see 
Sanders & McGaugh 2002 for a detailed review). In this case, the addition of dark 
matter to the "classical" mass is not required* and a typical spiral galaxy only consists 
of bare disk and bulge. 

Recently, Tiret & Combes (2007) presented a model of the Antennae galaxies, re- 
placing the Newtonian gravitation with MOND. The morphology obtained is in good 
agreement with observational data, as shown in Figure V.ll. 

The tidal field of galaxies within the framework of MOND is expected to be different 
than the one found with Newtonian dynamics. However, nothing excludes a priori the 
existence of compressive regions of similar shape or even intensity. To investigate this 
interesting point, a project in collaboration with H. Zhao (University of St Andrews) has 
been started, with the aims of (i) computing the tidal tensor thanks to MOND in the 
configuration of the Antennae obtained with Newtonian physics (e.g. those presented in 
Chapter II) and (ii) repeating the experiment on a fully MONDian model. Lot of tests 
need to be performed first and unfortunately, it is yet too soon to present the results 
here. 



h at galactic scale. Several studies show that DM is still needed in the MONDian prescription at the 
scale of galaxy clusters (see Sanders & McGaugh 2002). 
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Figure V.ll: Left: observations of the Antennae from Hibbard et al. (2001). Green 
represents stars and the gas is blue. Right: simulation with MOND by Tiret & Combes 
(2007). Stars are shown in yellow/red shades and the gas is still blue. (Adaptation of 
Figure 2 of Tiret & Combes 2007.) 
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Every new beginning comes from some other beginning 's end. 

— Seneca 



1 Conclusions 



During this study, we have focused on purely gravitational simulations of interacting 
galaxies and their tidal field. A collection of numerical tools has been developed to 
monitor the physical quantities related to this field. Thanks to it, we have shown that 
a special mode of the tides, called compressive mode, exhibits striking similarities with 
the properties of young substructures like star clusters and tidal dwarf galaxies. 

• The evolution of the mass in compressive mode is strongly correlated to the star 
formation rate. Peaks exist in both quantities when the interaction is at its climax. 

• The compressive regions are situated at the positions of observed clusters and tidal 
dwarf galaxies. 

• Their typical durations (~ 10 7 yr) are comparable to the age of these structures. 

• The corresponding energy (from ~ 10 46 to 10 57 erg) is sufficient to have a significant 
influence on the dynamical and chemical evolution of these objects. 

• The importance of the compressive mode decreases for mild, distant interactions, 
though it is of prime importance for mergers. 



It has also been shown that, if the details of the statistics of the compressive tides 
vary with the configuration of the galaxies and their encounter, these conclusions are 
robust over a wide range of parameters. 

It is now clear that the compressive tidal mode plays a non-negligible role in the 
formation and early evolution of star clusters. This demonstrates that a link exists 
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between the large scales (~ 10 1-2 kpc) where galactic dynamics set the morphology and 
kinematics of stellar and gaseous structures, and the smaller ones (~ 10 1-2 pc) for which 
the hydrodynamics takes over and determines the formation and chemical evolution of 
stellar associations. 

The details of such a link and the importance of the tidal field with respect to the 
other phenomena involved still need to be quantified thanks to high-resolution hydrody- 
namical simulations of giant molecular clouds but also young stellar associations. This 
next step requires a precise tuning of the parameters ruling the formation of stars and 
thus, will be done with extra-caution during the next months thanks to smoothed par- 
ticle hydrodynamics and adaptive mesh refinement runs. Furthermore, a study of the 
dynamical evolution of a cluster embedded in a tidal field is underway. We are now 
introducing initial mass function and stellar evolution to model a broad range of internal 
phenomena in addition to the external effect of the tides. These simulations will help 
us to pin down the survival or destruction conditions of a star cluster, during the first 
100 Myr of its life. 

Hopefully, all these projects will bring answers to many of the questions raised by the 
present work. 



2 Now and next 

In the continuation and aside from this work, I have conducted several projects related 
to the topic of this PhD. For some of them, the first results have just become available 
and are the starting points of very exciting explorations. 

2.1 Star formation in M51 

The Antennae galaxies have been chosen as observational reference because they are 
the closest major merger known. However, other mergers deserve as much attention as 
them, like for instance the grand design galaxy M51 and its companion NGC 5195. In 
this object, the tidal field has already played a crucial role in the formation of large scale 
structures like tails and spirals. It is then reasonable to question whether the smaller 
scales are also affected by the tides and in which fashion. 

To address these points, I use the smoothed particles hydrodynamics simulation pre- 
sented in Dobbs et al. (2010) to compute the tidal field and look for compressive regions. 
Figure 1 shows the very first results. Nodes of compressive modes exist along the main 
tidal tail, which suggests to push further the analysis of their properties, with respect 
to e.g. the Jeans mass or length. Interestingly, we also note that the spiral arms (which 
have been created during the encounter) are free of compressive zones. 
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Thanks to the presence of gas in the model, we will monitor the regions of star forma- 
tion and check for links with those of compressive tides. We will also study the evolution 
of the star formation rate through the major steps in the history of the interaction. 



2.2 iV-body simulation of the Stephan's Quintet 

What applies to two galaxies should work for five. To confirm this basic statement, I 
focused on Stephan's Quintet, a compact group of five galaxies that shows a very complex 
web of tidal structures and hydrodynamical features, e.g. the largest shock known in the 
Universe. Up to now, many scenarios based on observations in several wavebands have 
been proposed to explain each structure. In this project, I suggest a new approach that 
consists in splitting the interaction history into well-defined and separated events and to 
link them sequentially to get the global "storyboard" of the group. 

To verify this hypothesis and constrain the possible scenarios, I have performed more 
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Figure 1: Morphology of M51, 180 Myr after the first passage of the intruder. Blue 
represents the surface density of the stars and the gas. As usual, red shows the particles 
in compressive mode. The intruder NGC 5195 is a point whose mass is one third of the 
the one of M51. 
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than 3 000 simulations of the quintet to find a set of parameters that reproduces the mor- 
phology and the kinematics of the entire group (see Figure 2). This manual exploration 
of the parameter space has shown that several scenarios found in the literature are very 
unlikely while others allow to fairly reproduce the observations. These results will be 
published soon and hopefully will give new lights on the observations to better constrain 
their interpretation. 

The next step of this work will take these initial conditions and include the gaseous 
component needed to re-create the shock. This will be part of a larger project (gathering 
American and French observers and theoreticians) that has already started with the 
PhDs of Jeong-Sun Hwang in Iowa and Pierre Guillard in Paris (see Hwang et al. 2010; 
Guillard 2010). 

Obviously, the tidal field of the quintet will also be computed! 




Figure 2: iV-body simulation of the five members of Stephan's Quintet, plus the fore- 
ground galaxy in the lower part of the Figure. Tidal tails and inner structures match well 
the observations, despite some issues on the exact positions of the progenitors. (Adap- 
tation of Figure 10 of Renaud, Appleton, & Xu 2010.) 
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Overview ^ 

This Appendix briefly introduces some of the physical 
concepts or relations used in the main text. 



1 Density & potential 

When considering analytical definitions of mass profiles, it is extremely useful to get 
the gravitational potential from the density and vice versa. In this Section, we detail the 
derivations of the relations between the two quantities. 



1.1 From the density to the potential: Poisson's law 

The gravitational force on a particle of mass m at the position r = xe x + ye y + ze z 
due to an element of mass 5m situated at r' = x'e x + y'e y + z'e z is 

5F W = G m p-T^p 5m ' ( A - X ) 

which can also be written thanks to the density 8m = p(r') d 3 r'. By summing on all 
mass elements, one gets the total force 

F(r) = GmJ p(r') dV. (A.2) 

First, we notice that 

v ( 1 \ = ~ x)e x + (V' ~ y)&y + (/ - z )&z = r'-r 

| j» | J | j» 1 3 | j» 1 3 
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Therefore, 



F(r) = G m I V ( 1 ) p(r') dV 



r — r 



= m V U? y d 3 r'^j . (A.4) 

We now define the gravitational potential 

0(r) = —G [ d 3 r', (A.5) 
J \ r ' — r l 

and write 

F(r) = -m V0(r) (A.6) 

and 

V ■ F(r) = -m (V • V) <p(r) = -m V 2 0(r). (A.7) 

Back to Equation A. 2, we can take the divergence of the force and substitute it with 
the Laplacian of the potential 



V • F(r) = G m / V ■ ( — ) p(r') dV 



r' — r 



r' — r 



-V>(r) = G I V • ( , , |3 ) p(r') dV. (A.* 



Box A.l 



We note that 

„ , r' — r \ d ( x' — x 

V 



r' - r|V dx y[( x > - x )2 + ( y / _ y )2 + ( z / _ z )2]3/2 
a / y'-y 



+ 



+ 



dy \ [(x' - x) 2 + (y' - yf + (z' - z) 2 ^ 12 
d I z' — z 

9Z y[( X > - x )2 + (y> _ y)2 + ( z / _ ^)2]3/2^ 

3 3 _(x' - a?) + (y' - y) + (2/ - z) 



J j-./ j» 1 3 J j./ j» 5 

3 (V — r) ■ fr' — r) 

"I - 0- 



I . y I I f ^ j» 1 5 



Vr'/r. (A.9) 
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The integrand in Equation A. 8 is null for r' ^ r. Therefore, we can take the density 
out of the integral to get 



-V-o(r) = G i>{v) I V ( r ,'_ * ) d'r'. (A. 10) 



Then, we change the reference frame of the divergence from r to r': 



V 2 0(r) = G p(v) I V r , • ( ^— ^ ) dV. (A.ll) 



Thanks to the divergence theorem, we also have 



V-o(r): 6>(r) / X^T 3 d 2 s. (A.12) 



Here, we can write 

d?s= |r'-r|(r'-r)de, (A.13) 
where <iG is the infinitesimal solid angle. Then, we get 

V 2 0(r) = G p(r) j dQ, (A.14) 

and finally the Poisson's law 

V 2 0(r) = AttG p(r), (A.15) 

which becomes 



1 \d_ f 2 d<P(r) 
r 2 dr \ dr 

for spherically symmetric systems. 



4ttG p(r) (A.16) 



1.2 From the potential to the density 

To invert Poisson's law and derive the density from the gravitational potential, we 
consider spherically symmetric systems only. At a given radius r, we can consider two 
regions: a sphere (inside) and a infinite shell (outside). By superposition, we can write 

0(r) = in (r) + out (r) V r 6 [0, oo[. (A.17) 

Newton's first theorem stipulates that the force due to a shell of matter is null, everywhere 
inside the shell. Therefore, the potential there must be constant and we can estimate it 
at any point. E.g., <f> out {r) = ou t(O). 
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By replacing d 3 r' with 4nx 2 dx in the definition of the gravitational potential 
(Equation A. 5), and applying the formula for r = 0, one gets 



0out(V) = -47rG / xp out (x) dx. (A. 18) 

Jo 

However, by construction, p out (x < r) = 0. Thus, 

/oo 
xp ut(x) dx. (A. 19) 

Furthermore, Newton's second theorem tells us that the force from the inside region is 
equivalent to the one of all its matter concentrated in a point-mass at the center. We 
obtain this enclosed mass by summing the masses of infinitesimal shells 

M r) = -^<^> = f xV.nM dx. (A.20) 

r r J 

And finally, the global potential is 

0(r) = -4tcG (- J x 2 p(x) dx + J xp(x) dx^j . (A.21) 

Note that by applying the Laplace operator to this expression, and considering the limit 
0(oo) — > 0, one quickly finds Poisson's law. 



2 Dynamical friction 

In restricted simulations, a galaxy is represented by a point-mass surrounded by 
pseudo-particles. Therefore, the orbit of the galaxy follows Kepler's laws of motion. 
However, an interesting fact about galaxy mergers is that they actually merge! Both 
galaxies loose energy when they collide and spiral on each other until they finally merge. 
This loss of orbital energy is called dynamical friction. 

Consider an object flying through a population of stars uniformly distributed 
(Figure A.l). The closest stars are attracted by this object and tend to group around 
it. However, because of the proper motion of the flying object, this gathering of mass 
actually occurs behind it. This leads to a non-uniform distribution of the mass that drags 
the object backward, as a retarding wake. Finally, the motion of the object is slowed 
down: it looses kinetic energy (which is transfered to the nearby stars). 

For galaxy mergers, that implies that the motion of an intruder through a primary 
galaxy will be strongly affected by dynamical friction. The initial kinetic energy of the 
intruder will be lost and its orbit would bend. An intruder initially set on a Keplerian 
parabolic orbit (eccentricity e = 1) will effectively have a bounded orbit (e < 1) after it 
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•• •• *•• 

Figure A.l: An object (red) moves in a field of stars uniformly distributed (black). The 
closest stars (blue) are attracted toward the position of the object, but as it moves away 
(right), these stars group into a retarding wake and attract the object back, which slows 
it down. 

has experienced dynamical friction, through the material of another galaxy. Figure A. 2 
and Figure A. 3 compare the orbits of two Plummer (1911) spheres (blue) with the equiv- 
alent Keplerian orbit (red). The effect of dynamical friction in visible in the orbital decay 
of the Plummer spheres, showing that they have lost a fraction of their initial orbital 
kinetic energy (see also Figure A. 4). 

This phenomenon has been studied for many years. Chandrasekhar (1943) gave an 
mathematical expression for the velocity change the moving body undergoes: 

F dyn = -4n\n(A)^f? (erf(X) - ^e^ 2 ) v, (A.22) 

where v is the velocity of the body, M the mass considered (body + medium) and p the 
density of the surrounding medium. X = v/(ay/2) denotes the ratio of the velocity to 
the one of a Maxwellian distribution with dispersion a. The term In (A), called Coulomb 
logarithm (by analogy with the theory of plasmas), is a constant depending on the 
maximum impact parameter between the bodies and their typical relative velocity. Its 
usual value ranges from unity to ~ 10. 

Even if this mathematical description suffers from strong simplifications 3 it is possible 
to implement it in restricted simulations and thus correct the motion of the galaxies. 
Although, a proper description of the structures forming during a collision requires self- 
gravity of this regions. Therefore, instead of the fast restricted runs, consistent A^-body 
simulations are generally performed on powerful computers. In this case, each particle 
has a non-null mass and fully participates in the global potential, and the Keplerian 
solution is not valid anymore b . 



a For example, it neglects the self-gravity of the medium, in particular in the retarding wake that 
forms behind the moving body. In addition, the value In (A) is rather difficult to compute and is often 
arbitrary chosen. 

b An analytical solution exists for the two-body problem only, in the general case. 
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Figure A. 2: Orbit of the centers of mass of two Plummer spheres (N = 10 000, see their 
actual size in the thumbnail) initially set on a Keplerian orbit (of eccentricity e = 0.5, 
red) at the position marked with x . See Figure A. 3 for a zoom-in on the central region 
(dashed rectangle). 
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x 



Figure A. 3: Zoom-in on the central region of Figure A. 2. The pericenter passages are 
marked with •, + and *, respectively. Note that the two objects do not spiral around 
the same point at the end of the simulation. This is because they undergo a very close 
passage which scatters an important fraction of their particles, leading to an offset in the 
computation of their centers of mass. 



201 



Appendix A 



Basic concepts 



p 1 — 1 — 1 — 1 — i — 1 — 1 — 1 — 1 — i — 1 — 1 — 1 — 1 — i — 1 — 1 — 1 — 1 — i — 1 — 1 — 1 — 1 — i — 1 — 1 — 1 — <-^\ 
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500 1000 1500 2000 2500 3000 

time 

Figure A. 4: Evolution of the distance between the centers of mass of the Plummer 
spheres (see Figure A. 2). As in Figure A. 3, the symbols (•, +, *) correspond to the 
pericenter passages. The two configurations give out almost identical results up to the 
first pericenter passage. Afterward, the dynamical friction redistributes the orbital energy 
and the amplitude of the orbits decreases. The non-regular evolution of the distance after 
t = 2700 corresponds to the scatter of particles (mentioned earlier) after the very close 
passage at £ ~ 2650. 
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Overview 

This Appendix details the analytical derivations of the 
main dynamical quantities of several density profiles. 
The notations G, <ft, p, T, r, M and M(r) refer to the 
gravitational constant, the potential, the density, the 
tidal tensor, the radial coordinate, the total mass, and 
the mass enclosed in a sphere of radius r, respectively 



1 Classical cases 



For simplicity, all the core or characteristic radii are designed by Tq, while po represents 
the characteristic density Note however that the physical definition of these quantities 
may vary from one model to the other, as discussed in Section V.1.4 (page 173). 



1.1 Plummer 



We consider the Plummer (1911) potential given by 

GM 



Vq + r z 

Using Poisson's law (Equation A. 15, page 197), we get the density 



(b.i; 



2,Mr 2 n , . 

p(r) = a rrr^ - ( ) 

An (rj) + r 2 ) 1 
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To get the integrated mass, one can solve 



r 2 , , , 3M r 5 



M(r) = 4tt / xy(x) dx = — I dx (B.3) 



r Jo + 



(which is possible through the change of variable tan# = x/tq), or simply note that the 
system must obey Newton's second theorem 

= _ v « r) , (B.4) 

and get (much quicker!) 

Mr 3 
[Tq + r l ) 1 

We deduce that the radius r encloses ~ 35% of the total mass. 

Starting with the potential (Equation B.l) and using Einstein's summation conven- 
tion, we can derive the tidal tensor (see Equation 10, page 11): 

T " = GM9 ' 1 (B ' 6) 



and then 



T ij = -GM^ (rg + r2) -^ (B.7) 
(tq + r 2 ) ' 

(with = 1 if % — j and otherwise), which gives, for r = x (i.e. y = z = 0): 

T xx (r) = -GM—^ — . (B.8) 

(tq + r 2 ) ' 



1.2 Logarithmic potential 

This profile is defined by 

0(r) = ^ o 2 ln(l + £), (B.9) 

where vq denotes the circular velocity for the outer regions of the profile. The density is 
given by Poisson's law 

Pir) = ^ 2 f^l (B.10) 
(r 2 + rjj) 
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and the integrated mass comes from Newton's second theorem 

M(r) = v2 °gW^T) ■ (B11) 

We note that this mass integrates to infinity. Thus, it is impossible to use the term M 
in our calculations. For the tidal tensor, we have 

T- = -v^ (-J^) (B.12) 



and then 



which gives, for r = x 



,5 lj (r 2 + r 2 ,) — 2x { x 
(r 2 + rl) 2 



T» = ~v f " - " 2 > (B.13) 



^2 ^,2 



T xx (r) = vl- Sfa. (B.14) 



1.3 Hernquist 

The same can be done for the Hernquist (1990) profile: 

4>{r) = — • 

r + r 

One can easily find the density 



()'■' 



^"2 J iy 



(r + rf 

which gives, for r = x 

2GM 
(r + r) 

(see also Section V.2, page 175, in the framework of the 7-family with 7 = 1). 



(B.15) 



Mr 

P(r) = ~ ° , 3 - (B.16) 
lixr {ro + r) 

Once again, the integrated mass is found thanks to Newton's second theorem 

Mr 2 

M(r) = _ — (B.17) 
(r + r) 

We compute the tidal tensor 

T ij = -GMd i ( ^-2- J , (B.18) 

\(r + r) rj 

and get 



(B.19) 



T xx (r) = — ^ (B.20) 
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1.4 7-family 

Dehnen (1993) and Tremaine et al. (1994) denned the 7-family of density-potential 
pairs as follows: 



p(r) 



(3 - 7 )M r 

47T rT(r + r) 4 ' 1 



(B.21) 



and 



GM 



r (2- 7 ) 
GM 



r + r 



(2-7)' 



In 



r 



r + r 



if 7^2 



if 7 = 2, 



(B.22) 



which is equivalent to the Hernquist (1990) profile for 7 = 1, and to the Jaffe (1983) 
profile for 7 = 2. The integrated mass is therefore 

M r 3 " 7 



M(r) : ., . 

For 7 7^ 2, the associated tidal tensor is 

T l ^ 2 = -GMd i [ Xj r^ (r + r) 7 ~ 3 ] 



(B.23) 



GM 



r7 ( r o + r ) 



3-7 



rp . rf . ry* . <y» . 



(r + r) r 



(B.24) 



We note that the expression of Tj, 2 is equal to T^ =2 . So, we have the general expression 



GM 



3-7 



r7 ( r o + r ) 

for all 7's, which gives, along the x-axis 

T xx (r) = GM 



i:j _ 3r + r 7 
(r + rj r z 



2r + r (7 - 1) 



r7 ( r o + r 



,4-7 



(B.25) 



(B.26) 



(see also Section V.2, page 175). 



1.5 NFW 



The density of the NFW (Navarro, Frenk, & White 1997) profile is given by 



p(r) 



r 1 + 



2 • 



(B.27) 
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Using Equation A. 21 (page 198), we get 

/ 

0(r) = -4nGp r - 
\ 



x dx 

+ 



1+ x. 
ro 



dx 



\ 



\ r 



(B.28) 



/ 



Recall that the first integral (that we can solve with, e.g., a partial integration) is the 
integrated mass. Finaly, we obtain 



(r) = -47rGpor 2 -hi(l + - 



and 



M(r) = 47rp r 3 



In 1 + - 



r 7 r + r 



(B.29) 



(B.30) 



As for the logarithmic potential, the mass of an NFW profile integrates to infinity. The 
tidal tensor is given by 



T 3 = AirGfwl& 



x; In 1 + 



r 2 (r + r) 



(B.31) 



and finally 



T ij = A7rGp r 3 



(3xiXj r- r ° - 5 ij r" 3 ) In ^1 + -^J 



Ob %X j 



3 x ^ x j 



r 2 (r + r ) r 3 (r + r ) 2 ^ 4 (r + n>) 

which gives, for r = x 

T~(r) = (^) 
(see also Section V.2, page 175). 



21 , 1+ M 3r 2 + 2rr 



r / (r + r ) z 



(B.32) 



(B.33) 



1.6 Einasto 

Less classical than the previous ones, the Einasto profile (Einasto 1965; Cardone, 
Piedipalumbo, Sz Tortora 2005) has recently been presented as a better match to simu- 
lated dark matter halos (Navarro et al. 2010). The density-potential pair is defined by 



p(r) = p exp\ 

a 



r 



(B.34) 
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0(r) 



GM | r (a) ~ T («; Ij) /2 



- r (!) 



l/a 



ri'UrfV" 



a ar. 



o 



(B.35) 



where and T(x, y) are the gamma function and incomplete gamma function, defined 
by 



T(x) = r(x, 0) and r(x,j/)= / t x - L e- 1 dt. 

'y 



(B.36) 



The integrated mass is again derived from the potential thanks to Newton's second 
theorem 



M(r) = M 



with 



M 



a \ a I \a I \a 



[B.37) 



(B.38) 



Box B.l 



To compute the tidal tensor, we note that 
6 2r c 



<9 J r - 



a ar, 



-x^V-'/'VsV^exp 



2r c 



ark 



(B.39) 



Hence, 



GM 

rfW) 



<)' { —xf 3 



ro r(i) +Ij ^ a -/V-exp 



+XjT 6 r Q T -; 



3 2r a 



a ar. 



x^^a 1 ' 3 '^ 2 exp 



r (-) 1 



r 1 1 *1 



2r a 



a ar. 



-ii — 

ct 



Then, we have 



GM f 5 ij - 3r" 2 x,x 



r(i) 



3 2r c 



-XiXjr-^^a 1 - 3 ^ 3 exp 



2r a 



r 



2r a 



ark 



^B.40) 



(B.41) 
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Again, for r = x we have 



T xx (r) 



GM 



r 



3 2r a 



a ar, 



r 



- 2 3/a a 1 ~ 3/a rQ 3 exp 



2r c 



an 



(B.42) 



(see also Section V.2, page 175). 



2 Working with kernels 



When using iV-body simulations, a physical object is represented by a discrete dis- 
tribution of particles. To limit the effects of this numerical technique, the properties of 
each particle (in our cases, the gravitational potential) are smoothed over a scalelength 
e. The functional form that indicates the way of smoothing is the kernel. 

It is important to be familiar with the kernels when using them for the calculation of 
the gravitational force, or the tidal tensor: the kernel is explicitely used when substracting 
the effect of the particle in the center of a finite differences cube (see Section 1.2, page 38) 
and for the computation of the tensor itself in the alternative method (see Section 1.2.3, 
page 45). 

Dehnen's tree code (gyrf alcON, Dehnen 2000) proposes four different kernels, derived 
from density profiles. Kernel P0 is the classical Plummer form, detailled in the previous 
Section (with r = e). The next one, PI is a modified version of Plummer. Its density 
reads 



ppi(r) 



15M e 4 



8tt (e 2 + r 2 ) 



2^/2' 



(B.43) 



Using again Equation A.21, one gets 



15 GM 
2~T~ 



e r /e x 2 dx f°° x dx 

r./n (1 + X 2 ) ?/2 (1 -J- ^7/2 







A (l + x 2 



(B.44) 



The first integral can be solved by setting x = tan 9 and using trigonometry relations: 



r/e x 2 dx 
o (1 + x 2 ) 7 ' 2 



arctan(r/e) 



(cos 3 e - cos 5 e) do. 



(B.45) 
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Box B.2 



It is interesting to know that (V n > 0) 



J cos n 9 d9 = sin cos 71 - 1 + {n - 1) J sin 2 cos™" 2 d0 



cos"" 2 9 d9- I cos n 9 d9 



sin0 cos' 1 " 1 9 + (n - 1) 

sin cos™" 1 + (n - 1) J cos n " 2 rf0 



(B.46) 



We finally find 



r/e x 2 dx r 3 (5e 2 + 2r 2 ) 



/o (l + x 2 ) 7/2 15(e 2 + r 
Simpler, the second integral in Equation B.44 is 

x dx e 5 



2^/2 



which gives 



A (1 + x 



>pi{ r > 



2^/2 



5 (e 2 + r 



2^/2' 



GM (3e 2 + 2r 2 
2 (e 2 + r 2 ) 3/2 



From this simple form, it is easy to find the tidal tensor 
GM 



1 P1 



2 (e 2 + r 2 ) 7/2 



[x iXj (6r 2 + 21e 2 ) - 5 lj (2r 4 + 5e 4 + 7r 2 e 2 )] 



(B.47) 



(B.48) 



(B.49) 



(B.50) 



which can be implemented as r Jjf in the alternative method (see Equation 1.16, page 46). 
The potentials of the other kernels (not used in this work) can be derived similarly: 



pP2(r) 
4>P2{r) 



15M 
16tt 



7e 6 



2e 4 



(e 2 + r 2 ) 9/2 



e 2 + r 



2^/2 



GM (9e 4 + 10e 2 r 2 + 4r 4 ) 



4 (e 2 + r 



2^/2 



and 



pP3{r) 



105M 
64vr 



9e s 



4e 6 



(e 2 + r 2 ) 



2^1/2 



e 2 + r 



2^/2 



GM (45e 6 + 70e 4 r 2 + 56e 2 r 4 + 16r 6 



16 (e 2 + r 2 ) 



2^/2 



(B.51) 
(B.52) 

(B.53) 
(B.54) 
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Note that both P2 and P3 yield negative densities in their outer parts, corresponding 
to super- newtonian forces. Figure B.l displays the density profiles and the associated 
potentials. 




Figure B.l: Top: Density profiles (absolute value) of the four kernels implemented in 
gyrf alcON. The negative part (outer regions) of P2 and P3 is shown with dashed lines. 
Bottom: Associated potentials compared with the Newtonian solution (i.e. no smooth- 
ing). The same total mass has been used in all cases. 
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Appendix 

Glitches in the gravity 



Overview ^ 

In this Appendix, we investigate the sources of errors 
appearing in the numerical computation of the tidal 
tensor. 

) 



1 The problem 

In Chapter I, we have presented a finite differences method to compute the tidal ten- 
sor via the accelerations measured on a cube. Using a iV-body realization of a Plummer 
sphere, we have shown that the profile of the dominant term of the tensor is well re- 
produced, despite a certain level of noise. In the Figure 1.14 (page 44) however, some 
glitches are visible for well-defined radii. They are more visible in the zoomed-in version 
presented in Figure C.l. 

The errors mainly consist in a systematic, brutal shift from the expected correct 
solution to a lower value of the tidal term, over small radial ranges. Note that they do 
not exist in the alternative method (see Section 1.2.3, page 45) which simply sums the 
contributions of all individual particles to the tidal tensor. 

The origin of this problem has been investigated by several means, presented below. 



2 Possible sources 



The fact that other methods give solutions without these shifts indicates that the 
their origin is independent of the iV-body aspect of the system itself. Therefore, the 
rendition of the potential is exonerated from the list of possible causes. Remains the 
finite differences scheme and the computation of the acceleration. 
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Figure C.l: T xx term of the tidal tensor computed with the finite differences method 
and alternative method, compared with the analytical solution. 



The first order difference method implemented to compute the derivative of the ac- 
celeration introduces a numerical noise. However one would expect such a noise to be 
present everywhere in the computation and to yield a continuous evolution along, say, the 
radial axis. The phenomenon observed here is clearly different as it appears locally, and 
through discrete changes. That is, the glitches can not originate from the differencing 
scheme and therefore come from a problem due to the computation of the acceleration. 
Furthermore, the effect tends to disappear when the tolerance angle 9 goes to zero, as 
suggested on Figure C.2, and to completely disappear for a direct summation method 
(called alternative method, recall Figure C.l). 

Note however that a small tolerance parameter implies that more and more particles 
are considered for the direct summation computation, which weakens the power of the 
treecode. In the following (as well as for the rest of the study), we chose to sacrifice 
a small amount of precision in order to keep high computational performances, and set 



Hence, it is possible to narrow the field of possibilities to the tree aspect of the code 
and the way the accelerations are computed between the (sub) cells. To explore this 
in more details, one has to compare the value of the acceleration given by gyrf alcON 



e = o.4. 
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Figure C.2: Same as previous, but for different values of the tolerance parameter 9 used 
in the building of the tree. 



with those from other codes. In this matter, Figure C.3 displays the relative error in 
the evaluation of the acceleration, between Dehnen's code (gyrf alcON) and a classical 
treecode, hackcode. 

Errors are of the order of magnitude of 0.1% for the modulus of the acceleration 
(|a| = + + of). The same study on individual components (a x , a y or a z ) reveals 
that, for a very small number of particles (~ 10 over 10 5 ), these errors can reach 10 000%, 
but are compensated by a high accuracy for the other components. 

In the same time, Figure C.3 reveals a small, yet significant difference between the 
accelerations computed by gyrfalcON or getgravity (the tool used to compute the 
acceleration on the Active cube around each stellar particle). As parts of the falcON 
suite, both tools use the same libraries to create the tree and compute the accelerations, 
which makes the measured difference very strange. An explanation for this has not been 
found during this work and the problem had have to be overcome. 
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3 Partial solution 



Re- writing Dehnen's code was clearly not the goal of this study. Instead, it would have 
been possible to use the alternative method to compute the tidal tensor with the highest 
possible accuracy. Indeed, the only error introduced by this method is the softening kernel 
of the particles, a default which can not be avoided in iV-body simulations. However, 
the direct summation has a very high cost of CPU time and does not allow the statistical 
study (i.e. with a large number of particles and simulations) that we have conducted. 

The glitches in the getgravity method do not concern a large fraction of the parti- 
cles. In addition, it is likely (yet, impossible to check) that the problem does not exist 
for more complex configurations (where the tree is forced to be non-symmetric), like 
the ones used in this work. Finally, the glitches have a limited amplitude that may not 
compromise the determination of the sign of the eigenvalues, in other words, the com- 
pressive character of the tidal field. To be on the safe side however, we decided to set 
a "gray- zone" of width 0.14: if the maximum eigenvalue of the tidal tensor sits inside 
this range (—0.07 < A max < 0.07), the study of the tidal field is inconclusive. That is, 




Figure C.3: Relative error in the evaluation of the acceleration between gyrf alcON and 
hackcode (classical treecode). The solid lines indicates the median value of the relative 
error between the different tools. 
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the compressive tidal mode gathers the particles yielding a tidal tensor with a maximum 
eigenvalue less than -0.07 (in numerical units). 
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Overview 



This Appendix gives the detailed mathematical deriva- 
tion of the torque created by a galaxy on one of its 
clusters. Each Section presents a special configuration 
of both the galaxy and the cluster, starting from a sim- 
ple example and going to more general (and involved) 
cases. 



1 A dipole cluster around a point mass galaxy 

The goal of this derivation is to get a mathematical representation of the effect of a 
host galaxy on a non-spherical star cluster, say, a spheroid. A non-isotropic tidal field 
would then induce a torque on the cluster and possibly trigger (or enhance, or slow down) 
its rotation. 

To illustrate our approach, let's consider a textbook example where the cluster is 
made of two point-mass m\ = mi = m c /2 separated by a distance 2a along the vertical 
axis, and where the galaxy is a point mass mc situated at distances r% and r 2 from the 
cluster points (see Figure D.l). 

The potential of the cluster evaluated at the position of the galaxy is 



The distances are evaluated with trigonometric relations in spherical coordinates: 



J 






+ a 2 — 2arQ cos 9 
7q + a 2 + 2arQ cos 9, 



(D.2) 
(D.3) 
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Figure D.l: A cluster is represented by two points (blue) situated at distances r\ and r 2 
of a point mass galaxy G. Such a configuration has an axial symmetry, which makes it 
independent of (p, in spherical coordinates. 



and can be expanded in Legendre polynomials P n (see Equation D.ll below) to give 

i 1 00 / \ n 



i 1 / \ n 

- = — £ -- Pn(cos^) (D.5) 



^2 r G — ; V r G 



where r G > a, so that 



„ 00 / x 2n 

G m r \ / a 



^ „=o v u 



E - P 2n(cOS0). (D.6) 



By truncating the sum for n > 1, one gets the monopole and quadrupole terms (the 
dipole being zero): 



, , Gm c 
r G w 



(D.7) 



l + ^(3cos 2 0- 1 
2r c 



The torque induced by the galaxy on the cluster about the origin (i.e. the center of mass 
of the cluster) is therefore 

r G / c = -r c/G w -m G r G x [-V0 c (r G )] 

Gm c mc a 2 



cos sin (D.8) 



Note that for a = 0, one retreives the case of a monopole cluster which yields no torque. 

Such a result can be obtained because of the simplicity of the setting. This corresponds 
well to the case of globular clusters which exhibit a regular shape. For younger objects 
like proto-clusters which have more filamentary structure, a more detailed derivation is 
required and the use of spherical harmonics is preferable. 
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2 A cylindrical cluster 

We now write the expression of the torque for a mass distribution of the cluster 
symmetric with respect to the 9 = axis (i.e. depending on the radial and latitudinal 
but not longitudinal coordinates). 



Box D.l 



We adopt here the Condon-Shortley convention, as in Binney & Tremaine 
'1987). A spherical harmonic is defined by: 



*tW>V>)-y 47r {t+\ m \)\ F t {COSV)e X \l form<0 ' 

(D.9) 

where Pj m ' is the associated Legendre function of the first kind defined by: 

TO^i-T"^ (D.10) 

with P £ ° the Legendre polynomial, itself defined by: 

1 d e (x 2 -l) e 

m-) = Pe^) = ~ (X dx , l> - (D.11) 
For negative degrees, one has the relation 

*r B (^v) = (-i) m *r(^v). (D-i2) 

where the asterisk denotes the complex conjugation of the harmonic ty. 



Consider a spheroidal cluster of potential c (r) embedded in a galaxy represented by 
a density pc( r )- The net torque of the galaxy on the cluster is the opposite of the one 
from the cluster on the galaxy, given by the time derivative of the angular momentum 
L G : 



Tg/c = — r c / G 



dL, 



dt 



[ p G (r)[rx V<f) c (r)]dV. (D.13) 
Jv 



In our approach, both the density of the galaxy and the potential of the cluster are 
decomposed in spherical harmonics on concentric spherical shells, to ease the evaluation 
of their gradients. Anderson (1992, see also Kawai & Makino 1998) proposed a fast 
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multipole method (or FMM) to express the potential of the cluster at the position r, 
from its potential c (os) on the surface of a sphere of radius a < r: 

1 00 +1 f 

71 n=0 r J S r 

where P n is the n-th Legendre polynomial, and ds = sm8 s d9 s dip s . We decompose c (as) 
in spherical harmonics: 



, . Gm c 
las) — 



a 

k=0 



J2C k *Ue s ,p s ), (D.15) 



(The upper index of \I/ translates the independence of the potential on ip s .) The 
argument of the Legendre polynomial in Equation D.14 is (s • r)/r = cos a, with a 
the angle between direction vectors s and r which define the (9 s ,(p s ) and (9, ip) angles, 
respectively . This can be re-writen as 

A/jr n 

P n (cosa) = —j *nW><P) (D.16) 



2n + 

j=-n 



Therefore, one has 



's ■ r\ , , N , Gm c 4tt 



J P n <fic(as)ds 



a 2n + l 



00 n „ 

E E °* / *2(*-,<P.) da. (D.17) 



fc=0 j = — n 

The orthogonality between the spherical harmonics implies 



/ y° k (e s , v% ) ¥:{e B , Vt ) ds = (D.18) 

Js 



and thus 

E 

n=0 

Computing the gradient, one has 



^^E0"^*i(^)- (D.19) 



rxV «r) = -^i + M) # (D . 20) 
sin at/? at' 

= — 2jW C " 90 v ' (D ' 21) 

n=0 
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and 



86 



2n + 1 dP n (cost 



d6 



2n + 1 „ dPJ cos 6) 
smfc> 



47T 



dcosd 



2n + 1 

47T 



Pi f COS I 



(D.22) 



The index = 1 of the associated Legendre function implies that the term n = vanishes. 
Physically, this means that the monopole term of the potential of the cluster does not 
participate in the torque induced by the galaxy. Finally, one gets 



Gm r 



G/c 



E 



71=1 



Pg(tJ 



a\ n+1 



P„(cos9) <p dV. 



(D.23) 



2.1 First case: point mass galaxy 



The simplest case is to consider the host galaxy as a point-mass at the position r G . 
Thus, the density reads pc( r ) — m G <H r ~~ r c)- Equation D.23 becomes 

n+l 



G/c 



Gm c m G 



E 



n=l 



2n + l o / _o_ 

4tt " Vr G 



(D.24) 



Back to our textbook example (Section 1), we can identify the potential of the cluster 
at r G to its decomposition in spherical harmonics (Equation D.19) up to the quadrupole, 
with the help of Box D.l: 



rc V 47r 



Arc 



Tq V 16-7T 



(D.25) 



and get, by comparison with Equation D.7, 



°0 



C° = 




(D.26) 



By using these values in Equation D.24, one retrieves the result of Equation D.8. 



2.2 Second two point galaxy 



If now the galaxy is made of two points symetrically distributed with respect to the 
cluster, i.e. at the positions (r G , G , <p G ) an d (r G , 7r — 9 G , tt + <p G ) the density becomes 



m G 5(r - r G ) 



71 



e. 



— — — 6{<p - <p G ) + — - 

sin 9 G sin i- Uc 



Gj 



6{(p-n- ip G ) 



. (D.27) 
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Such a configuration has new symetrical properties, which allow us to verify our calcu- 
lations. 

■ Box D.2 



One has to be careful when using Dirac functions in spherical coordinates. 
In this system, a volume element is not rectangular because of the Jacobian 
matrix used in for the transformation of coordinates from the cartesian base. 
This implies that the point r G = (r G , 9 G , ip G ) is not represented by 

S(t - r G ) ^ 8(r - r G )5(6 - 9 G )5(<p - <p G ), (D.28) 

but by 

6(r - r G ) 5(9 - 6 G 
q sin 9 G 

so that the 3D Dirac is normalized to unity 



5(r - r G ) = v 2 w v . w %> - ^g), (D.29) 

Tn Sin 9t 



[ 8(r-r G )dV = [ 5(r - r G ) r 2 dr sin 6(19 dip = 1. (D.30) 
Jv Jv 



With this density, Equation D.23 becomes 

r 



OA- = -^pi± ^/^C» (± Y +I {P„'[co S » G ] - P„'[oo S (n - @g)] J <Pq- 

n=l ^ ' 

(D.31) 

One immediately has -P^[cos (it — 9 G )] = P^[— cos (9 G )] an d, knowing that 

Pi{-x) = (-lR 1 ^), (D.32) 

one can see that the even terms cancel out because of the symmetry of the galaxy. Finally, 
one gets 

T G/C = _^£f J^±C< L+1 (^) 2 " +2 pi, +1 (oo S « G ) a . (D.33) 



2.3 General case for the galaxy 

It is also possible to represent the density of the galaxy with a new set of spherical 
harmonics, evaluated on a shell, without loss of generality: 

p G (r) = '" G %~ ' C) E £ W«4 (D.34) 



G 



£=0 m=-£ 
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so that Equation D.23 becomes 



Gm c m G \- \- \- 2n + 1 

n=l 1=0 m=-l 

y (-) 5(r-r G )vI/-(^^)P n 1 (cos0) £ dV. (D.35) 

One can write 

*?(0,¥>) = K™P l e ml (cos 0) e im ^, (D.36) 

with 



21+1 (£-jmj)l f form>0 
A£ "V 4vr (£+| m |)! X \l form<0 ' 

and thus 

(-) 5(r -re)® 1 ? {8, <p) P^{cos 6) <p dV =K™ j (-J <5(r - r G )r 2 

P^(cos0)PJ m| (cos0)sin0 d0 







/>2tt 

/ e %mtp (p dip. (D.38) 



The integral on r reads 



/•°° , a \ n+l { n \ n+1 

J (7) 6 ( r - dr = r c [yj > ( D - 39 ) 



and the one on 9 can be re-writen 



P r ](cos 6) P 1 ™ 1 (cos 6) sin 6 d6 = J P r ](u)P l e ml (u) du. (D.40) 
The last integral is best expressed in the (^-independent cartesian base: 

/*27T /•27T 

/ [cos (mip) + 1 sin (mip)\ [cos yy — sin yx] dip = y / cos (my) cos y <iy 

/>2tt 

— x / cos (my) sin y ciy? 

/>2vr 

+ ty I sin (my) cos ip dip 
Jo 



2- 







— zx / sin (my) siny dy. (D.41) 
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We note that only the first and the last integrals are non-zero, giving 

/■27T 

/ [cos (rrup) + i sin (imp)] [cos ipy — sin (pit] dip = n [(y + zx) <5 m 1 + (y — zx) 5^] . 
•/ii 

(D.42) 

Therefore, |m| = 1 which forbids £ = 0, so that the monopole term of the galaxy does not 
create a torque on the cluster, in the reference frame centered on the cluster 3 . Further- 
more, the orthogonality condition on the associated Legendre functions in Equation D.40 
gives 

= ^l< n + (D.43) 

Back to Equation D.35, one obtains 

rt+l 



n=l m=— n 



iC^T [(y + **) C + (y - «*) <SJ . (dm) 

Furthermore, we recall from Equation D.37 that 



K- 1 = -K\ = , 1 1V (D.45) 

4tt n(n+l) 



To derive the coefficient we note that 



mG ^ 2 ~ rG) ^. (D.46) 



By applying this relation backward to D n q *, one gets 



= (-iy ™ G %~ rG) Di 

r G 

D - r = (_i)ff£,«. (D.47) 



as oposed to the usual reference taken at the galactic barycenter. 
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In particular, one has 

D- 1 = -D 1 :, (D.48) 

so that 

D- l K~ l = D%Kl (D.49) 

Thus, one finally has 

r-t 00 / \ n+1 

Gm c m G 



G/c 



a 

n 



f \ n+ 

E V^+l) - C° [9l(Di)y + , (D.50) 



where and denote the real and imaginary parts of the complex z. It is also 

possible to define an ang le $ so that tan/3 = -3(Z£)/9t(l£), and write 



g/c = ^£^ £ v^TT) f-V" C° (cos/3 y + sin/3 x) 

^^^f v^TT) c: \Dl\ ft, (D.51; 



where \D\\ is the modulus of the complex number D\. 



3 General case 



Up to now, we have considered the cluster to be cylindrically symmetric. If now we 
go to the general case, the gradient takes a more complicated form. By analogy with 
Equation D.19, the potential of the cluster evaluated on the sphere of radius a < r is 
now 



r =- 



, EE©" <%nv,<p)- (D.52) 

n=0 q=—n 

We still have 



sm.6 dip o9 

but now, both terms are non-zero: 



a# a 

n=0 q=—n 



EE(") ™ t ^ (D .54) 

n=0 q=—n 



eta a 

n=0 q=—n 
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One has 

dPn ( cos 1 

d6~ 



— [ sin" 



sin 191 (0) 



S q P%{cos6 
dcosM (6) 



d9 

g sin 191 1 (6) cos - n ; . ,J + sin 1 " 1 (6) 
11 v 7 dcosl'l (0) K ' 

f*OS f7 i 

\q\—EPl q \™sO)-P\«\ + \ C os9). 
sin u 



S q P2[ 



COS( 



d9 dcos6\ dcosM (0) 

(D.56) 



Note that this expression is not defined for q = ±n because of the second associated 
Legendre function. However, for these values, one can demonstrate that 



dPn (COS I 

d8~ 



i ,COS0 

9 — ~E P n cos( 
sin 6* 



(D.57) 



i.e. same as previous but without the second term. Finally, 



Gm r 



G/c 



Mr 



E E c ;(; 



n=0 <?= — n 



r 



iq 



sin 6 1 



**(0,^) 0+|? 



COS T Gf / /l \ a 

^ *«.(«.¥>>¥> 



n-1 



- E 

g=— n+l 



a\ n + 1 Kl 



K 



_ ^(9-191-1)^1+1(0^) ^ I dK 



(D.58) 



From Equation D.37, one recalls that 
Kl 



kl- 



- = v^+ 



n — q l 



l?l x 



■1 for g > 

_l)kl+i forg<0 



(D.59) 



and one retrieves the result of Equation D.23 for q = (i.e. cylindrical symmetry). 

It is now possible to replace the density of the galaxy with its decomposition in 
spherical harmonics (as in Section 2.3). Though technically possible, such an analytical 
derivation requires many terms to give a precise description of a real galaxy. Therefore, 
the use of iV-body simulations would be prefered to a purely analytical work. In this 
case, the volume integral can be replaced by a discrete sum over all the particles of the 
galaxy: 



Gm c 



N 



G/c 



E^EIE^ 

Kl 



j=l n=0 \q=—n 

n-1 , ^ n+l 

E 



n+l 



1(1 nvs^Oi+M^nViMWs 

olll Uj 



sin 6j 



(D.60) 



q=— n+l x 

where rrij is the mass of the j-th particle of the galaxy, located at coordinates (r 3 -, 9j, ipj). 
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4 Opposite torque 

The derivations presented above use the potential of the cluster and the density of 
the galaxy to compute the torque. In some cases, it is easier to consider the opposite. 
One could be tempted to simply switch the two subscripts and get minus the torque 
(r c /G = —Tg/ c ). If this is true in principle, one should not forget that the way the 
potential is evaluated from the Anderson's sphere depends on whether one stands inside 
or outside this sphere. Equation D.14 gives 



for r > a, which is true when the potential of the cluster is needed at the position of 
the galaxy. But when doing the opposite, the potential of the galaxy is required at the 
positions of the stars of the cluster, i.e. inside the sphere (r < a). In this case, one has 
to write 



Because of the cross product with r, the gradient of 4> along r is never computed and 
thus, this change does not affect the final results more than transforming 



i 00 +i r 

n=0 JS 



(D.61) 



1 p 

^) = ^ E( 2n + X ) (~) / P « (v") ^(«s) ds 



(D.62) 



n=0 




into 




(D.63) 



when switching the subscripts G and c. 
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Appendix E 

Virial theorem 



Overview 

In this Appendix, we derive the virial theorem in the 
general case, and then in the context of an isotropic 
tidal field. We also present the concept of virial radius. 



1 Derivation of the scalar virial theorem 

1.1 General case 

Let's consider a system of N particles. Its moment of inertia is given by 

N 

/ = $> 4 rf. (E.l) 



By deriving it, we obtain 



dl 
~dt 

2dt 2 



N 



2 £ 

i=i 

N 

^rrii di n + m i v i 
i=i i=\ 

N N 

Y,F l r t + Y,^vl (E.2) 



i=i 

N N 



i=l i=l 



where is the acceleration of the i-th particle. The forces and their potential energies 
can be writen as 

E = k and F — -V 'E = -k £ (E.3) 



233 



Appendix E 



VlRIAL THEOREM 



where k is a constant and £ an integer. By replacing the term (Fj 7™j) in Equation E.2, 
one gets 

potentials 

where is the kinetic energy. This equation is called general virial theorem but is rarely 
used. In the case of a stationnary system, the first derivative of the moment of inertia is 
null, and a fortiori is the second. Therefore we have a simpler form for the theorem: 

2K- £E = 0. (E.5) 

potentials 



1.2 Gravitational case 

Usual assuptions in stellar dynamics set that the only interaction in the system is the 
gravitation (£ = —1), whose potential energy is Q. Therefore, the virial theorem for a 
stationary system (Equation E.5) with gravitational interaction only becomes 

2K + n = 0. (E.6) 



1.3 In a tidal field 

Let's consider a star cluster embedded in an isotropic tidal field represented by the 
tensor 

/ A \ 

T = [0 A . (E.7) 
\ A / 

The differencial force between an element % of mass and the barycenter of the cluster 
is 

5F tidefM = rmX rj. (E.8) 
We note that in this case £ = 2 and thus 

^tides,; = -^ m i^ r i- ( E - 9 ) 

When summing on all the mass elements, if we assume that the tidal field is uniform 
over the entire cluster, we have 

-Ftides,i n = / Ar 2 dm 

= Xa MR 2 t , (E.10) 
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where M is the total mass of the cluster, R t is the radius where the density drops to 
zero, and 

i r M 

a =W„ r dm (E - n) 

(see Section E.3 for various examples). Finally, we can write for the entire cluster 

£tides = -\\a M Rl (E.12) 
and, from Equation E.10, the new virial theorem becomes 

2K + Q + Xa M R 2 t = 0. (E.13) 



2 Velocity dispersion and virial radius 

In the relation above, global quantities have been derived by summing over all the 
mass elements of the cluster. It is also possible to use quantities that globaly characterize 
the cluster. For example, the velocity dispersion a can be used in the expression of the 
kinetic energy: 

1 - 1 
K = -Y j m l vl = -Ma' (E.14) 

and the same can be done with the gravitational energy 

« = (E,5) 

where r v is called virial radius a . Therefore, the total energy of a system in virial equilib- 
rium (without tidal field: K = — fl/2) is 

E = K + Sl= V = _GM 1 (E . 16) 
2 4r v K J 



3 a parameter 

As defined in Equation E.ll, the a parameter describes the mass distribution of a 
spatially extended object. Using usual mass profiles (like, e.g. those of Appendix B), it 
is possible to calculate a via the density p and the relation dm = Anxp{r)r 2 dr: 

47T f°° 

a= MRl.L Pir)rdr - (EU7) 

a This definition is based on Meylan & Heggie (1997), but Binney & Tremaine (1987) set another 
expression. 
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Note that the upper limit of the integral can now be set to infinity because the density 
will drop to zero where the integrated mass reaches M. 



3.1 Homogeneous sphere 

Let's first consider a homogeneous sphere of density p and radius i? t 



Immediately, one gets 



and, with 



one has 



« = (E.19) 

5M * v ' 



Po = ^L, (E.20) 



a = % (E.21) 
5 



3.2 Power law profile 



Let's now consider a density profile defined by 



We have 



p(r) = { P° fe) 7 if r < ^ (E.22) 
else. 



Mi£ 2 



t J o 

7 d3-7 



M(5-7) 



We also know that 



M = Anpor^ I r 2 ~ 7 dr 
Jo 



3-7 



so, 



5 — 7 

In the case 7 = 0, we retrieve the homogeneous sphere. 



(E.23) 



(E.24) 



a = (E.25) 
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3.3 Plummer sphere 

For a Plummer sphere, we use the density given in Equation B.2 (page 205) and 
change the variable thanks to tan# = r/ro, to get 



a 



3r 2 f°° r 4 dr 



Rl Jo (r 2 + r 2 ) 5/2 
3r 2 f$ tan 4 9 dO 



R t Jo C os 2 #(^) 

V cos^ 1 



5/2 



3r 2 f e sin 4 9d9 



R 2 J cos 9 
3r 2 f 6 sin 4 fl 



R( Jo 1 - sin' 9 



cos 9 d6. (E.26) 



Then, we set y = sin 9 and get 



a 3r 2 [V y* dy 



Jo l-l/ 2 

?2 



r> >2 ^argtanhy - - - yj . (E.27) 

However, if the Plummer model yields a finite mass, it does not have a truncation ra- 
dius. Therefore, we have to set an artificial one. Knowing the integrated mass (see 
Equation B.5, page 206), we can see that 99% of the total mass M is enclosed in the 
sphere of radius Rt ~ 12r . 



Box E.l 



To determine the value of y: 

y = sin 9 = sin 



arctan I — 
r 



(E.28) 



let's use the complex notation (i 2 = —1 and 3t(z) is the imaginary part of the 
complex z). We know that 



sin z = 3 [exp {iz)\ and arctan z— — In f — — — 

2i \ 1 — iz 



(E.29) 
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Therefore, 



sin (arctanz) = £y I W — - — 1 = —=^= (E.30) 



1 - %Z ) y/T+ 

which is obviously true for real numbers too. 



With our arbitrary truncation radius, we get 

y=-^L= (E.31) 

which gives 

a « 3.86 x 10" 2 . (E.32) 

A similar development can be used for the Hernquist profile, this time by setting the 
arbitrary truncation radius at R t = 200r . In this case, one obtains a = 9.33 x 10 -3 . 

4 Star formation efficiency 

4.1 In isolation 

The total energy of young cluster in isolation is simply the sum of its potential and 
kinetic energies: 

E = K + Q 

This cluster formed with a SFE e defined as the ratio between the stellar mass to the 
total (stars + gas) mass: 

Now, let's assume that this cluster expels its gaseous content, via stellar feedback. If 
the expulsion is rapid enough (impulsive mass loss, see also Hills 1980; Boily & Kroupa 
2003a), one can consider that neither the velocity dispersion nor the virial radius have 
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changed during the process. Therefore, the energy becomes 




(E.35) 



Assuming that the cluster was initially virialized, one has 



Ma 2 = 



GM 2 
2r v ' 



(E.36) 



therefore 



Ei 



GM 2 e 

4r v 
GM 2 e 

2r v 




GM 2 e 2 
2r v 



(E.37) 



The stars of such a cluster remain bound if the energy is negative, i.e. if e > 50%. In 
other words, when a molecular cloud collapses, its material is gravitationally bound in a 
potential well created by a mass M. When the first stars blow the remaining gas away, 
this mass changes, which flattens the well and allows the less bound stars to escape. If 
the fraction of the remaining gas is less than 1/2, the mass loss is not important enough 
to allows this unbinding process, and the cluster survives 



When the star cluster is not isolated but embedded in a tidal field, one must take the 
tidal energy into account: 



Impulsive approximation 

Once again, after the instantaneous expulsion of the gas (which leaves a, r v and Rt 
unchanged), one can write 



4.2 In a tidal field 



Eq = K + Vt + -Etides 




V 



(E.38) 




GM 2 



-\a M* R 2 t 



2r v 




(E.39) 
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Figure E.l: Schematic view of the energy profile of a cluster. With no tidal field (green), 
the criterion for a bound system is E < 0. When adding an extensive field (blue), this 
level is shifted so that even a cluster with E < can be unbound. This situation does 
not occurs with a compressive tidal field (red). 



Using the virial theorem (Equation E.13) at the initial stage, one gets 

E x = Ge M2 (I - e \ _ e \ a M R 2 . (E.40) 



2r v V 2 

Because of the tidal term, the potential energy is now the sum of the intrinsic potential 
of the cluster, and the external potential (via the tides) seen as an harmonic potential 
(Figure E.l). This shifts the reference level of energy and thus, a bound system does not 
necessarily correspond to E < 0. Instead, we compare the energy E\ of the cluster after 
the gas expulsion to those of a (virtual) cluster of the same mass but virialized, and still 
embedded in the same tidal field: 

E' = -KUa 2 -\a' M* R' 2 

2 2r' v 2 * 

= e -Mo 2 - ^J^l _ i Xa ' M R'l (E.41) 
2 2r' v 2 *' v ; 

which becomes, thanks to the virial theorem, 

E' = - \ -e Xa' M R'l (E.42) 

K 

Hence, our cluster remains bound if its energy after expulsion E\ is lower than the 
reference level E' of a virialiazed system of same mass, i.e. 

Ei < E' (E.43) 



240 



4 - Star formation efficiency 



G e M 2 
4r v 



1 - 2e + e — - e \ M [a R; - a R A < 



(E.44) 



We set 



A = e Xa M R 2 



GeM 2 \~ l XaR 2 r v 



GM 



(E.45) 



which leads to 



1 - 2e + e ( 



AA 



a \R t 



< 0. 



(E.46) 



Adiabatic evolution 



If one does not consider the expulsion of gas as an instantaneous process but rather 
as a slow evolution, the energy of the gas-free cluster takes into account the modified 
shape parameters (a', r' v and -R' t ): 



Ey = -M^a 2 l X a> M* R'l. 

2 2r, 2 * 



(E.47) 



Again, we use the virial theorem from the initial state: 



Ex 



G e M 2 ( 1 



2r v 



-e A M ( a R 2 t + a' R' 2 t ) . 



(E.48) 



This energy has to be compared to the enegy of the same system if it was virialized 
(recall Equation E.42): 



G e M 2 
4r v 



(l - e ^ - X -e A M (a R 2 t - a' R'fj < 0. 



(E.49) 



With Equation E.45, we obtain 



1-e 



2A 



1 - 



a' [R' t 



a \ R 



< 0. 



(E.50) 



Solution 



In both cases (impulsive or adiabatic), one is able to get the properties (mainly the 
radii) of the bound cluster after the gas expulsion, depending on the strenght of the 
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tidal field (which is encapsulated in the term A). However, the analytical solving of 
Equation E.46 or Equation E.50 is only feasible with a simplifying assumption of the 
ratios a' /a, r' v /r v and R' t /R t . 

The minimal hypothesis which consists in assuming a homologous transformation (i.e. 
r' v /r v = R' t /Rt) does not recover the different responses of the core and the external 
part of the cluster to the mass- loss, as noted by Boily & Kroupa (2003b). However, 
an analytical solution does not exist in the general case, because the relations giving 
the evolution of the radii and mass profile are not known. An example of derivation in 
the frame of the (non-realistic) homologous transformation is given in the Section III. 2. 2 
(page 116). 
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Overview 

This Appendix lists the publications (in preparation, 
submitted, in press, refereed and conference proceed- 
ings) I collaborated to. A short teaser (instead of an 
abstract) is available for the major papers. 



1 In prep., submitted... 



• Renaud F., Boily CM. & Theis C. 

Dark matter halos profiles and tidal fields 
MNRAS Letters, submitted 

In this short contribution (5 pages, 3 figures), we study analytically the tidal field of 
DM halos profiles, some derived from cosmological simulations (NFW, Einasto). Our 
results show that the superposition of some potentials (Dehnen's 7-family, Einasto), 
as in configurations close to galaxy-galaxy interactions, yields large areas of compres- 
sive tides, while NFW do not. Because the star formation is very likely enhanced 
when embedded in a compressive tidal field, we claim that the precise shape of 
the DM halo should be chosen with care, when modeling star formation in galaxy 
mergers. 

• Renaud F., Appleton P. & Xu K. 

N-body simulation of the Stephan's Quintet 
ApJ, submitted 

In this paper (12 pages, 13 figures), we propose a possible scenario for the formation 
of the Stephan's Quintet (a compact group of five galaxies) thanks to collisionless 
iV-body simulations. With an innovative serial approach, we cover a large range 
of parameters and rule out several scenarios. We discuss the possible future of 
the group, in the more general frame of the hierarchical formation and evolution 
scenarios for compact groups of galaxies. 
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2 Refereed papers 

[1] Karl S., Naab T., Johansson P., Kotarba H., Boily CM., Renaud F. & Theis C. 

One moment in time - Modeling star formation in the Antennae 
ApJ Letters, 715, L88 [ADS if ] 

In this Letter (6 pages, 4 figures), we present a SPH simulation of the Antennae 
galaxies. The hydrodynamical aspects of star formation, feedback and radiative 
cooling are included to measure the SFR along time, and the distributions of star 
clusters in the central region. An enhancement of the SFR is noted shortly after the 
two passages. This new model complements our purely gravitational approach by 
confirming the role of the large-scale interactions on the demographics of the star 
clusters. 

[2] Renaud F., Boily CM., Naab T. & Theis C 
Fully compressive tides in galaxy mergers 
2009, ApJ, 706, 67 [ADS & } 

This paper (16 pages, 16 figures) is the direct continuation of the Letter Renaud et 
al. (2008) published one year before in MNRAS. The concept of compressive tides is 
introduced analytically and explored through a series of iV-body simulations. The 
Antennae galaxies are taken as an observational reference. Then, a survey on the 
parameters of interacting galaxies is presented, to broaden the conclusions from the 
Antennae. The paper closes with a discussion on the link between compressive tidal 
field and star formation, in term of energy, enhancement of the SFR and multiple 
stellar populations. 

[3] Renaud F., Boily CM., Fleck J.-J., Naab T. & Theis C 

Star cluster survival and compressive tides in Antennae-like mergers 
2008, MNRAS Letters, 391, L98 [ADS * ] 

In this Letter (5 pages, 3 figures), we briefly introduce the concept of compressive 
tides and concentrate on an iV-body model of the Antennae galaxies. Details of the 
spatial distribution of the compressive regions, as well as their characteristic times 
are presented and confronted to recent observational data on young clusters in this 
merger. 
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3 Conference proceedings 

[4] Renaud F., Theis C, Naab T. & Boily CM. 

Substructures formation induced by gravitational tides? 
Galaxy Wars - Johnson City Tn, USA (19-22 July 2009) 
2010, ASP Conference Series, 423, 191 [ADS if ] 

[5] Hwang J.-S., Struck C, Renaud F. & Appleton P. 

Models of the intergalactic gas in Stephan's Quintet 
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